max-norm optimization problem (Definition 3.2). Then there exists a quantum
algorithm that outputs an $\tilde{x}\in[0,1]^{n}$ satisfying (3.16) with
$\epsilon=1/3$ using $O(\sqrt{n}\log n)$ quantum queries to $O_{f}$, with
success probability at least $0.9$.

In other words, the quantum query complexity of the max-norm optimization
problem is $\tilde{\Theta}(\sqrt{n})$.

We prove Theorem C.1 also using search with wildcards (Theorem 3.1).

###### Proof.

It suffices to prove that one query to the wildcard query model $O_{c}$ in
(3.1) can be simulated by one query to $O_{f_{c}}$, where the $c$ in (3.14) is
the string $c$ in the wildcard query model.

Assume that we query $(T,y)$ using the wildcard query model. Then we query
$O_{f_{c}}(x^{(T,y)})$ where for all $i\in[n]$,

$\displaystyle x_{i}^{(T,y)}=\begin{cases}\frac{1}{2}&\text{if }i\notin T;\\\
0&\text{if }i\in T\text{ and }y_{i}=0;\\\ 1&\text{if }i\in T\text{ and
}y_{i}=1.\end{cases}$ (C.38)

If $c_{|T}=y$, then

  * •

if $|T|=n$ (i.e., $T=[n]$), then

$\displaystyle f_{c}(x)=\max_{i\in[n]}|x_{i}^{(T,y)}-c_{i}|=0$ (C.39)

because for any $i\in[n]$, $x_{i}^{(T,y)}=y_{i}=c_{i}$;

  * •

if $|T|\leq n-1$, then

$\displaystyle
f_{c}(x)=\max_{i\in[n]}|x_{i}^{(T,y)}-c_{i}|+g_{i}=\frac{1}{2},$ (C.40)

because for all $i\in T$ we have $x_{i}^{(T,y)}=y_{i}=c_{i}$ and hence
$|x_{i}^{(T,y)}-c_{i}|=0$, and for all $i\notin T$ we have
$|x_{i}^{(T,y)}-c_{i}|=|\frac{1}{2}-c_{i}|=\frac{1}{2}$.

Therefore, if $c_{|T}=y$, then we must have
$f_{c}(x^{(T,y)})\in\big{\\{}0,\frac{1}{2}\big{\\}}$.

On the other hand, if $c_{|T}\neq y$, then there exists an $i_{0}\in T$ such
that $c_{i_{0}}\neq y_{i_{0}}$. This implies $x_{i_{0}}^{(T,y)}=1-c_{i_{0}}$;
as a result, $f_{c}(x^{(T,y)})=1$ because on the one hand
$f_{c}(x^{(T,y)})\geq|1-c_{i_{0}}-c_{i_{0}}|=1$, and on the other hand
$f_{c}(x^{(T,y)})\leq 1$ as $|x_{i}^{(T,y)}-c_{i}|\leq 1$ for all $i\in[n]$.

Notice that the sets $\big{\\{}0,\frac{1}{2}\big{\\}}$ and $\\{1\\}$ do not
intersect. Therefore, after we query $O_{f_{c}}(x^{(T,y)})$ and obtain the
output, we can tell $Q_{s}(T,y)=1$ in (3.1) if
$O_{f_{c}}(x^{(T,y)})\in\big{\\{}0,\frac{1}{2}\big{\\}}$, and output
$Q_{s}(T,y)=0$ if $O_{f_{c}}(x^{(T,y)})=1$. In all, this gives a simulation of
one query to the wildcard query model $O_{c}$ by one query to $O_{f_{c}}$.

As a result of Theorem 3.1, there is a quantum algorithm that outputs the $c$
in (3.14) using $O(\sqrt{n}\log n)$ quantum queries to $O_{f}$. If we take
$\tilde{x}=c$, then $f_{c}(\tilde{x})=\max_{i}|c_{i}-c_{i}|=0$, which is
actually the optimal solution with $\epsilon=0$ in (3.16). This establishes
Theorem C.1. ∎


# Jupiter’s Dynamical Love Number

Dong Lai Cornell Center for Astrophysics and Planetary Science, Department of
Astronomy, Cornell University, Ithaca, NY 14853

###### Abstract

Recent observations by the Juno spacecraft have revealed that the tidal Love
number $k_{2}$ of Jupiter is $4\%$ lower than the hydrostatic value. We
present a simple calculation of the dynamical Love number of Jupiter that
explains the observed “anomaly”. The Love number is usually dominated by the
response of the (rotation-modified) f-modes of the planet. Our method also
allows for efficient computation of high-order dynamical Love numbers. While
the inertial-mode contributions to the Love numbers are negligible, a
sufficiently strong stratification in a large region of the planet’s interior
would induce significant g-mode responses and influence the measured Love
numbers.

dynamical tides — giant planets — Jupiter’s interior — gravitational fields

##  1 Introduction

The Juno spacecraft recently found an “anomaly” in Jupiter’s tidal Love
number: the measured $k_{2}=0.565\pm 0.006$ (Durante et al. 2020) appears to
be smaller than the theoretical hydrostatic value $k_{2}^{\rm(hs)}=0.590$
(Wahl et al. 2020) by $4\%$. This discrepancy may be explained in terms of
dynamical tides, i.e., Jupiter’s response to the finite-frequency tidal
forcings from the Galilean moons (Idini & Stevenson 2021). Here we present a
simple calculation that explains this Love number “anomaly” quantitatively.
Naive expectation would suggest a $1/(\omega_{\alpha}^{2}-\omega^{2})$
enhancement (where $\omega_{\alpha}$ is the f-mode frequency of the planet) of
the tidal response due to the finite tidal frequency ($\omega$) as compared to
the hydrostatic ($\omega=0$) response. The key to obtain the correct answer is
to treat the rotational (Coriolis) effect on the modes of a rotating planet
and their tidal responses in a self-consistent way. Our general method also
allows for efficient computation of high-order dynamical Love numbers
$k_{lm}$, as well as the inclusion of the contributions to $k_{lm}$ from the
inertial modes (due to planetary rotation) and g-modes (due to stable
stratification in the planetary interior).

##  2 Dynamical Love Number and Normal Modes

Consider a planet (mass $M$, radius $R$ and spin angular frequency
$\Omega_{s}$) orbited by a satellite (mass $M^{\prime}$) in a circular orbit
with semi-major axis $a$ and orbital frequency $\Omega_{\rm orb}$. We assume
the spin axis is aligned with the orbital axis. In the frame corotating with
the planet, the $(lm)$-component of the tidal potential produced by
$M^{\prime}$ on the planet is

$U({\bf r},t)=-A_{lm}\,r^{l}\,Y_{lm}(\theta,\phi)\,e^{-i\omega t},$ (1)

where $A_{lm}=(GM^{\prime}/a^{l+1})W_{lm}$ (with $W_{lm}$ a dimensionless
constant; $W_{lm}\neq 0$ when $l+m=$even), ${\bf r}=(r,\theta,\phi)$ specifies
the position vector (in spherical coordinates) measured from the center of the
planet, and

$\omega=m(\Omega_{\rm orb}-\Omega_{s})$ (2)

is the tidal forcing frequency. It suffices to consider only $m>0$. The
relevant non-zero tidal components are $(lm)=(2,2),(3,1),(3,3),(4,2),(4,4)$
etc.

The linear response of the planet to the tidal forcing is specified by the
Lagrangian displacement, ${\mbox{\boldmath$\xi$}}({\bf r},t)$, of a fluid
element from its unperturbed position. In the rotating frame of the planet,
the equation of motion takes the form

$\frac{\partial^{2}{\mbox{\boldmath$\xi$}}}{\partial
t^{2}}+2{\bf\Omega}_{s}\times\frac{\partial{\mbox{\boldmath$\xi$}}}{\partial
t}+{{\bf C}}\cdot{\mbox{\boldmath$\xi$}}=-\nabla U,$ (3)

where ${\bf C}$ is a self-adjoint operator (a function of the pressure and
gravity) acting on $\xi$ (see, e.g., Friedman & Schutz 1978). A free mode of
frequency $\omega_{\alpha}$ (in the rotating frame) with
${\mbox{\boldmath$\xi$}}_{\alpha}({\bf
r},t)={\mbox{\boldmath$\xi$}}_{\alpha}({\bf r})\,e^{-i\omega_{\alpha}t}\propto
e^{im\phi-i\omega_{\alpha}t}$ satisfies

$-\omega_{\alpha}^{2}{\mbox{\boldmath$\xi$}}_{\alpha}-2i\omega_{\alpha}{\bf\Omega}_{s}\times{\mbox{\boldmath$\xi$}}_{\alpha}+{\bf
C}\cdot{\mbox{\boldmath$\xi$}}_{\alpha}=0,$ (4)

where $\\{\alpha\\}$ denotes the mode index, which includes the azimuthal
number $m$. We carry out phase-space mode expansion (Schenk et al. 2002)

$\left[\begin{array}[]{c}{\mbox{\boldmath$\xi$}}\\\
{\partial{\mbox{\boldmath$\xi$}}/\partial
t}\end{array}\right]=\sum_{\alpha}c_{\alpha}(t)\left[\begin{array}[]{c}{\mbox{\boldmath$\xi$}}_{\alpha}({\bf
r})\\\ -i\omega_{\alpha}{\mbox{\boldmath$\xi$}}_{\alpha}({\bf
r})\end{array}\right].$ (5)

Using the orthogonality relation
$\langle{\mbox{\boldmath$\xi$}}_{\alpha},2i{\bf\Omega}_{s}\times{\mbox{\boldmath$\xi$}}_{\alpha^{\prime}}\rangle+(\omega_{\alpha}+\omega_{\alpha^{\prime}})\langle{\mbox{\boldmath$\xi$}}_{\alpha},{\mbox{\boldmath$\xi$}}_{\alpha^{\prime}}\rangle=0$
(for $\alpha\neq\alpha^{\prime}$), where $\langle
A,B\rangle\equiv\int\\!d^{3}x\,\rho\,(A^{\ast}\cdot B)$, we find (Lai & Wu
2005)

${\dot{c}}_{\alpha}+i\omega_{\alpha}c_{\alpha}={iQ_{\alpha,lm}\over
2\varepsilon_{\alpha}}\,A_{lm}\,e^{-i\omega t},$ (6)

where

$\displaystyle
Q_{\alpha,lm}\equiv\bigl{\langle}{\mbox{\boldmath$\xi$}}_{\alpha},\nabla(r^{l}Y_{lm})\bigr{\rangle}=\int\\!d^{3}x\,r^{l}\,Y_{lm}\,\delta\rho^{\ast}_{\alpha},$
(7)
$\displaystyle\varepsilon_{\alpha}\equiv\omega_{\alpha}+\langle{\mbox{\boldmath$\xi$}}_{\alpha},i{\bf\Omega}_{s}\times{\mbox{\boldmath$\xi$}}_{\alpha}\rangle,$
(8)

and we have used the normalization
$\langle{\mbox{\boldmath$\xi$}}_{\alpha},{\mbox{\boldmath$\xi$}}_{\alpha}\rangle=1$.
In Eq. (7), $\delta\rho_{\alpha}$ is the Eulerian density perturbation
associated with the eigenfunction ${\mbox{\boldmath$\xi$}}_{\alpha}$.

Equation (6) has stationary solution

$c_{\alpha}(t)={Q_{\alpha,lm}\over
2\varepsilon_{\alpha}(\omega_{\alpha}-\omega)}\,A_{lm}\,e^{-i\omega t}.$ (9)

The gravitational perturbation associated with the density perturbation
$\delta\rho({\bf r},t)=\sum_{\alpha}c_{\alpha}(t)\delta\rho_{\alpha}({\bf
r})$, evaluated at the planet’s surface ($r=R$), is

$\delta\Phi({\bf r},t)\Bigr{|}_{r=R}=-\sum_{\alpha}c_{\alpha}(t)\,{4\pi\over
2l+1}\,{GMQ_{\alpha,lm}\over R}\,Y_{lm}.$ (10)

Thus the tidal Love number is

$k_{lm}={\delta\Phi\over U}\Bigr{|}_{r=R}={2\pi\over
2l+1}\sum_{\alpha}{\bar{Q}_{\alpha,lm}^{2}\over\bar{\varepsilon}_{\alpha}(\bar{\omega}_{\alpha}-\bar{\omega})}.$
(11)

In the above equation, the tidal overlap coefficient $\bar{Q}_{\alpha,lm}$ and
the mode frequencies $\bar{\omega}_{\alpha}$ and $\bar{\varepsilon}_{\alpha}$
are in units where $G=M=R=1$, i.e.,
$\bar{\omega}_{\alpha}=\omega_{\alpha}/(GM/R^{3})^{1/2}$, etc. Note that for a
given $m>0$, the sum in Eq. (11) includes modes with positive
$\omega_{\alpha}$ and negative $\omega_{\alpha}$, corresponding to prograde
(with respect to the planet’s rotation) and retrograde modes.

##  3 F-mode Contribution

Table 1: Oscillation modes of non-rotating polytropic $(n=1)$ planet model

|  | $\omega_{0}$ | $Q_{l}$ | $C$
---|---|---|---|---
$\Gamma_{1}=2$ |  |  |  |
$l=2$ | f | 0.1227E+01 | 0.5579E+00 | 0.4991E+00
| p1 | 0.3462E+01 | 0.2690E-01 | 0.1119E+00
$l=3$ | f | 0.1698E+01 | 0.5846E+00 | 0.3321E+00
| p1 | 0.3975E+01 | 0.4054E-01 | 0.8627E-01
$l=4$ | f | 0.2037E+01 | 0.5979E+00 | 0.2489E+00
| p1 | 0.4409E+01 | 0.4623E-01 | 0.6984E-01
$\Gamma_{1}=2.4$ |  |  |  |
$l=2$ | f | 0.1230E+01 | 0.5580E+00 | 0.4985E+00
| g1 | 0.4688E+00 | -0.1313E-01 | 0.1057E+00
| g2 | 0.3270E+00 | 0.3071E-02 | 0.1342E+00
| g3 | 0.2526E+00 | -0.8961E-03 | 0.1464E+00
$l=3$ | f | 0.1703E+01 | 0.5850E+00 | 0.3317E+00
| g1 | 0.5681E+00 | -0.1216E-01 | 0.3321E-01
| g2 | 0.4117E+00 | 0.3257E-02 | 0.5622E-01
| g3 | 0.3254E+00 | -0.1031E-02 | 0.6603E-01
$\Gamma_{1}=2.4$ | for | $r\in[0.5,0.7],$ | $[0.85,0.93]\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ |
$l=2$ | f | 0.1228E+01 | 0.5580E+00 | 0.4988E+00
| g1 | 0.3664E+00 | -0.6986E-02 | 0.9279E-01
| g2 | 0.2022E+00 | 0.1408E-02 | 0.1207E+00
| g3 | 0.1565E+00 | -0.1546E-03 | 0.1467E+00
$l=3$ | f | 0.1701E+01 | 0.5849E+00 | 0.3319E+00
| g1 | 0.4489E+00 | -0.6067E-02 | 0.2790E-01
| g2 | 0.2749E+00 | 0.2049E-02 | 0.4166E-01
| g3 | 0.2126E+00 | -0.1567E-03 | 0.6622E-01

Note. — $\omega_{0}$ and $Q_{l}$ are the mode frequency and tidal overlap
coefficient (Eq. 7), both in units such that $G=M=R=1$, and $C$ is defined in
Eq. (14). The planet’s density profile is that of $n=1$ polytrope (with the
equation of state $P\propto\rho^{2}$). The first model has $\Gamma_{1}$ (the
adiabatic index) equal to $\Gamma=1+1/n$, and we list the properties for the
f-mode and the first radial-order p-mode. The second model has
$\Gamma_{1}=2.4$ throughout the planet, and the third model has
$\Gamma_{1}=2.4$ only in two regions ($r/R\in[0.5,0.7],\,[0.85,0.93]$) and
$\Gamma_{1}=\Gamma$ otherwise (the transition width is $0.025R$; see Eq. 26),
and we list the properties for the f-mode and the first three radial-order
g-modes. Note that when $|Q_{l}|\ll 1$, the quoted $Q_{l}$ values are only
accurate in 2-3 significant figures.

In most situations, the sum in Eq. (11) is dominated by f-modes since they
have the largest tidal overlap $Q_{\alpha,lm}$. For planetary rotation rate
$\Omega_{s}$ much less than the breakup rate $(GM/R^{3})^{1/2}$, (e.g.,
$\bar{\Omega}_{s}=0.288$ for Jupiter), the effect of rotation on the modes can
be treated perturbatively (e.g. Unno et al. 1989). Let $\omega_{0}$ ($>0$) be
the mode frequency of a non-rotating planet, then for a given $m>0$, the sum
in Eq. (11) includes

$\displaystyle\varepsilon_{\alpha}\simeq\pm\omega_{0},$ (12)
$\displaystyle\omega_{\alpha}=\pm\omega_{0}-mC\Omega_{s},$ (13)

with

$\displaystyle mC$ $\displaystyle\equiv$
$\displaystyle\int\\!d^{3}x\,\rho\,{\mbox{\boldmath$\xi$}}^{\ast}_{\alpha,0}\cdot(i{\bf\hat{z}}\times{\mbox{\boldmath$\xi$}}_{\alpha,0})$
(14) $\displaystyle=$ $\displaystyle m\int_{0}^{R}\\!\\!dr\,\rho
r^{2}(2\xi_{r}\xi_{\perp}+\xi_{\perp}^{2}),$

where
${\mbox{\boldmath$\xi$}}_{\alpha,0}=\left[\xi_{r}(r){\bf\hat{r}}+\xi_{\perp}(r)r\nabla\\!_{\perp}\right]Y_{lm}$
is the mode eigenvector of a non-rotating planet. To a good approximation, we
can also set $Q_{\alpha,lm}$ to be the non-rotating value, i.e,

$Q_{\alpha,lm}\simeq Q_{l}.$ (15)

Thus Eq. (11) reduces to

$k_{lm}\simeq\left({4\pi\over
2l+1}\right){\bar{Q}_{l}^{2}\over\bar{\omega}_{0}^{2}-(mC\bar{\Omega}_{s}+\bar{\omega})^{2}}.$
(16)

Figure 1:  Dynamical correction $\Delta
k_{lm}/k_{lm}^{\rm(hs)}=(k_{lm}-k_{lm}^{\rm(hs)})/k_{lm}^{\rm(hs)}$ to
Jupiter’s tidal Love number as a funtion of the orbital frequency $\Omega_{\rm
orb}$ of the perturbing satellite (in units of the spin frequency
$\Omega_{s}$). All results (solid curves) are computed using the $n=1$
(isentropic) polytrope model, except that the dot-dashed curve is for
$k_{2}=k_{22}$ computed using the $n=0.9$ polytrope model. The vertical dashed
lines specify the orbital frequencies of Io, Europa, Ganymede and Callisto
(from right to left). The $-4\%$ “anomaly” of $k_{2}$ observed by Juno can be
explained by the planetary model with $n\simeq 1$.

For an incompressible planet model ($n=0$ polytrope), the $l=2$ mode (Kelvin
mode) has

$\bar{\omega}_{0}={2\over\sqrt{5}},\quad\bar{Q}_{2}=\left({3\over
2\pi}\right)^{1/2},\quad C={1\over 2}.$ (17)

Thus

$k_{2}\equiv k_{22}\simeq{3\over 2}\,\Bigl{[}1-{5\over
4}(\bar{\Omega}_{s}+\bar{\omega})^{2}\Bigr{]}^{-1},$ (18)

with $\bar{\omega}=2(\bar{\Omega}_{\rm orb}-\bar{\Omega}_{s})$.

Giant planets are approximately described by a $n=1$ polytrope (corresponding
to $P\propto\rho^{2}$). Table 1 list the numerical values of
$\omega_{0},\leavevmode\nobreak\ Q_{l}$ and $C$ for several non-rotating
poltropic models (with different levels of stratification; see Section 5). For
$l=m=2$ tidal response (and $n=1$), $2C\simeq 1$, we have

$k_{2}\simeq
0.520\,\Bigl{[}1-0.664(\bar{\Omega}_{s}+\bar{\omega})^{2}\Bigr{]}^{-1}.$ (19)

Applying to the Jupiter-Io system: Jupiter has $\bar{\Omega}_{s}=0.288$ [with
spin period 9.925 hrs and $2\pi(R^{3}/GM)^{1/2}=2.863$ hrs], Io has
$\bar{\Omega}_{\rm orb}=0.0674$ (orbital period 1.769 days), so
$\bar{\omega}=-0.441$. Thus the hydrostatic and dynamical $k_{s}$ values are

$k_{2}^{\rm(hs)}=0.550,\quad k_{2}=0.528=0.960\,k_{2}^{\rm(hs)}$ (20)

This explains the $4\%$ discrepancy between $k_{2}$ and $k_{2}^{(\rm hs)}$.
Note that our static $k_{2}^{(\rm hs)}$ does not agree with the value (0.590)
from Wahl et al. (2020). This could arise for two reasons: (i) The simple
$n=1$ polytropic model does not precisely represent Jupiter’s internal
structure; (ii) In deriving Eq. (16), we have neglected order
$\bar{\Omega}_{s}^{2}$ corrections to the mode frequency and the tidal
overlap111Both $\bar{\omega}_{0}^{2}$ and $\bar{Q}_{l}^{2}$ in Eq. (16) can
have corrections of order $\bar{\Omega}_{s}^{2}$. These corrections do not
affect $\Delta k_{2}/k_{2}^{\rm(hs)}$ to the leading order. Note that for
$n=0$ (incompressible MacLaurin spheroid), the exact expressions for the mode
frequency and tidal overlap coefficient are available [see Eq. (3.4) and Eq.
(3.22) of Ho & Lai (1999)]..

Figure 1 shows the dynamical corrections $\delta_{lm}\equiv\Delta
k_{lm}/k_{lm}^{\rm(hs)}=(k_{lm}-k_{lm}^{\rm(hs)})/k_{lm}^{\rm(hs)}$ to
Jupiter’s tidal Love numbers as a function of the orbital frequency
$\Omega_{\rm orb}$ of the perturbing satellite. These results are obtained
using the $n=1$ isentropic model, which gives the hydrostatic values
$k_{lm}^{\rm(hs)}\simeq 0.550,\,0.213,\,0.219,\,0.121,\,0.123$ for
$(lm)=(22),\,(31),\,(33),\,(42),\,(44)$, respectively. Although these
hydrostatic values may not correspond to the “true” values for Jupiter because
of the simplicity of the polytrope model and the $\bar{\Omega}_{s}^{2}$
corrections (see above), the dynamical corrections $\delta_{lm}$ shown in Fig.
1 are robust.

Our results depicted in Fig. 1 can be compared to those of Idini & Stevenson
(2021) obtained using more complicated calculations (see their Table 2). Our
$\delta_{22}$, $\delta_{33}$ and $\delta_{44}$ values (evaluated for the
orbital frequencies of Io, Europa, Ganymede and Callisto) agree reasonably
well with theirs, but our $\delta_{31}$, $\delta_{42}$ values are a factor of
a few smaller.

Finally, using Table 1, we can easily check that the contributions from
p-modes to $k_{\rm lm}$ are negligible.

##  4 Inertial-Mode Contribution

In addition to f-modes and p-modes, a rotating planet possesses a spectrum of
inertial modes supported by Coriolis force.

For $n=1$ polytrope, the $m=2$ inertial modes have been computed by Xu & Lai
(2017) using a spectral code. The mode properties are

$\omega_{+}=0.556\Omega_{s},\leavevmode\nobreak\
\varepsilon_{+}=0.28\Omega_{s},\leavevmode\nobreak\
\bar{Q}_{+}=0.015\bar{\Omega}_{s}^{2}$ (21)

for the prograde mode, and

$\omega_{-}=-1.10\Omega_{s},\leavevmode\nobreak\
\varepsilon_{-}=-0.55\Omega_{s},\leavevmode\nobreak\
\bar{Q}_{-}=0.010\bar{\Omega}_{s}^{2}$ (22)

for the retrograde mode, where $Q_{\pm}$ is the tidal coupling coefficient
$Q_{\alpha,22}$. Since $\omega<0$, we can write the inertial mode contribution
of $k_{2}$ as

$k_{2,\rm in}={2\pi\over
5}\left[{\bar{Q}_{+}^{2}\over\bar{\varepsilon}_{+}(\bar{\omega}_{+}+|\bar{\omega}|)}+{\bar{Q}_{-}^{2}\over|\bar{\varepsilon}_{-}|(|\bar{\omega}_{-}|-|\bar{\omega}|)}\right].$
(23)

Define $\hat{\omega}\equiv\omega/\Omega_{s}$ (and similarly
$\hat{\omega}_{\pm}$ and $\hat{\varepsilon}_{\pm}$) and
$\hat{Q}_{\pm}\equiv\bar{Q}_{\pm}/\bar{\Omega}_{s}^{2}$, we have

$k_{2,\rm in}={2\pi\over
5}\,\bar{\Omega}_{s}^{2}\left[{\hat{Q}_{+}^{2}\over\hat{\varepsilon}_{+}(\hat{\omega}_{+}+|\hat{\omega}|)}+{\hat{Q}_{-}^{2}\over|\hat{\varepsilon}_{-}|(|\hat{\omega}_{-}|-|\hat{\omega}|)}\right].$
(24)

For $n=1$ polytrope, this gives

$k_{2,\rm in}={2\pi\over 5}\,\bar{\Omega}_{s}^{2}\times
10^{-4}\left({8.04\over 0.556+|\hat{\omega}|}+{1.82\over
1.10-|\hat{\omega}|}\right).$ (25)

For Jupiter-Io system, $\hat{\omega}=2\Omega_{\rm orb}/\Omega_{s}-2=-1.53$, it
is clear that $k_{2,\rm in}\ll 1$. In general, unless $|\hat{\omega}|$ happens
to be very close to $|\hat{\omega}_{-}|$ (to within $10^{-4}$), the
contribution of the inertial modes to the Love number is negligible.

##  5 Stable Stratification and G-Mode Contribution

Figure 2:  Same as Fig. 1, but for the $n=1$, $\Gamma_{1}=2.4$ planetary
model, which possesses g-modes. The heavy solid curves include the
contributions of f-modes and first three radial-order g-modes to $k_{lm}$, the
light solid curves include only f-modes (for the $n=1$ isentropic model, as in
Fig. 1). Figure 3:  Same as Fig. 2, but for the $n=1$ model with stable
stratification only in the region $r/R\in[0.5,0.7]$ (see Eq. 26). Figure 4:
Same as Fig. 2, but for the $n=1$ model with stable stratification in two
regions: $r/R\in[0.5,0.7]$ and $[0.85,0.93]$.

In Sections 3-4 we considered fully isentropic models for Jupiter, i.e., the
adiabatic index $\Gamma_{1}\equiv(\partial\ln P/\partial\ln\rho)_{s}$ equals
the polytropic index $\Gamma\equiv d\ln P/d\ln\rho=1+1/n$. In reality, some
regions of the planet may be stably stratified, with $\Gamma_{1}>\Gamma$.
Indeed, the gravity measurement by Juno and structural modeling suggest that
Jupiter have a diluted core and a total heavy-element mass of 10-24 Earth
masses, with the heavy elements distributed within an extended region covering
nearly half of Jupiter’s radius (Wahl et al. 2017; Debras & Chaberier 2019;
Stevenson 2020). The composition gradient outside the diluted core would
provide stable stratification, and the planet would then possess g-modes.
Another stable region may exist between (0.8-0.9)$R$ and $0.93R$ (Debras &
Chaberier 2019).

To explore of how g-modes influence the tidal love numbers, we consider three
simple planetary models, all having a $n=1$ density profile ($\Gamma=2$), but
with different adiabatic index profiles: (i) $\Gamma_{1}=2.4$ throughout the
planet (see Table 1); (ii) $\Gamma_{1}=2.4$ only in the stable region
$\bar{r}=r/R\in[0.5,0.7]$ (with a transition width of 0.025) and
$\Gamma_{1}=\Gamma=1+1/n$ otherwise, i.e.,

$\Gamma_{1}(r)=2+{0.4\over[1+e^{40(\bar{r}-0.7)}][1+e^{-40(\bar{r}-0.5)}]}$
(26)

(iii) $\Gamma_{1}=2.4$ only in two stable regions $\bar{r}=r/R\in[0.5,0.7]$
and $[0.85,0.93]$ (with a transition width of 0.025) and
$\Gamma_{1}=\Gamma=1+1/n$ otherwise (see Table 1).

For each model, we compute the f-modes and g-modes of a nonrotating planet
(see Table 1), and use Eqs. (12)-(13) to account for the effect of rotation on
the modes. We include only the first three radial-order g-modes in our
calculation of $k_{lm}$. The perturbative approach of the rotational effect is
approximately valid for these modes since $mC\Omega_{s}$ is less than the mode
frequency $|\omega_{0}|$.

Figures 2-4 show the results for the dynamical Love numbers based on the three
models. It is obvious that significant dynamical correction to the hydrostatic
$k_{lm}^{\rm(hs)}$ occurs around the resonance, where
$\omega_{\alpha}=\omega$. The “strength” of each resonance is measured by the
tidal overlap coefficient, and a large $|Q_{l}|$ value implies that the
“width” of the resonant feature is larger (see Eq. 16). For the
$\Gamma_{1}=2.4$ model, the stratification is strong, the broad/strong
resonance with the g1 mode can affect $k_{lm}$ associated with the Galilean
moons (Fig. 2). For the model with the stable stratified region restricted to
$r/R\in[0.5,0.7]$ (Fig. 3), the resonance feature is much weaker/narrower, but
still $\Delta k_{2}/k_{2}^{\rm(hs)}$ becomes $-5\%$ for Io. When the model
further includes the stratified region at $r/R\in[0.85,0.93]$ (Fig. 4), the
resonance features shift and broaden, and $k_{31}$ becomes affected for the
Galilean moons. Obviously, these results are for illustrative purpose, but
they indicate that resonance features due to stable stratification in the
planet’s interior may influence the the measured dynamical Love numbers.

Note that Figures 2-4 do not include contributions from high-order g-modes.
These modes (with mode frequencies comparable to $\Omega_{s}$) become mixed
with inertial modes (so-called “inertial-gravity” modes; see Xu & Lai 2017)
and cannot be treated using Eqs. (12)-(13). However, because of their small
tidal overlap coefficients, they are unlikely to be important contributors to
$k_{lm}$ except for the coincidence of an extremely close resonance.

##  6 Conclusion

We have derived a general equation (Eq. 11) for computing the dynamical Love
number $k_{lm}$ of a rotating giant planet in response to the tidal forcings
from its satellites. In most situations, the Love number is dominated by the
tidal response of f-modes, and the general expression reduces to Eq. (16),
which can be easily evaluated using the mode properties of nonrotating planet
models (see Table 1). We show that the $4\%$ discrepancy between the measured
$k_{2}$ of Jupiter and the theoretical hydrostatic value can be naturally
explained by the dynamical response of Jupiter’s f-modes to the tidal forcing
from Io – the key is to include the rotational (Coriolis) effect in the tidal
response in a self-consistent way. We also show that the contributions of the
inertial modes to the Love number $k_{2}$ are negligible.

We have also explored the effect of stable stratification in Jupiter’s
interior on the Love numbers. If sufficiently strong stratification exists in
a large region of the planet’s interior, g-mode resonances may influence the
dynamical Love numbers associated with the tidal forcing from the Galilean
moons. Thus, precise measurements of various $k_{lm}$ could provide
constraints on the planet’s interior stratification.

## References

  * Bate (2009) Debras, F., & Chabrier, G. 2019, Astrophys. J. 872, 100
  * Bate (2009) Durante, D., et al. 2020, Geophysical Research Letters, 47, e2019GL086572
  * Bate (2009) Friedman, J.L., & Schutz, B.F. 1978, Astrophys. J. 221, 937
  * Bate (2009) Ho, W.C.G., & Lai, D. 1999, MNRAS, 308, 153
  * Bate (2009) Idini, B., & Stevenson, D.J. 2021, PSJ, in press (arXiv:2102.09072)
  * Bate (2009) Lai D., & Wu, Y. 2006, Phys. Rev. D 74, 024007
  * Bate (2009) Schenk, A.K. et al. 2002, Phys. Rev. D 65, 024001
  * Bate (2009) Stevenson, D.J. 2020, AREP, 48, 465
  * Bate (2009) Unno, W., et al. 1989, Nonradial oscillations of stars, University of Tokyo Press
  * Bate (2009) Wahl, S. M., et al. 2017, Geophys. Res. Lett. 44, 4649
  * Bate (2009) Wahl, S. M., Parisi, M., Folkner, W. M., Hubbard, W. B., & Militzer, B. 2020, The Astrophysical Journal, 891, 42
  * Bate (2009) Xu, Wenrui, & Lai, D. 2017, Phys. Rev. D 96, 083005


Snehal Rajput, Rupal Agravat, Mohendra Roy, Mehul S. Raval 11institutetext:
Department of Information and Communication Technology,
Pandit Deendayal Petroleum University, Gandhinagar, India.
22institutetext: School of Engineering and Applied Science,
Ahmedabad University, Ahmedabad, India.
22email<EMAIL_ADDRESS>
22email<EMAIL_ADDRESS>22email<EMAIL_ADDRESS>
# Glioblastoma Multiforme Patient Survival Prediction

Snehal Rajput 11 Rupal Agravat 22 Mohendra Roy 11 Mehul S. Raval 22

###### Abstract

Glioblastoma Multiforme is a very aggressive type of brain tumor. Due to
spatial and temporal intra-tissue inhomogeneity, location and the extent of
the cancer tissue, it is difficult to detect and dissect the tumor regions. In
this paper, we propose survival prognosis models using four regressors
operating on handcrafted image-based and radiomics features. We hypothesize
that the radiomics shape features have the highest correlation with survival
prediction. The proposed approaches were assessed on the Brain Tumor
Segmentation (BraTS-2020) challenge dataset. The highest accuracy of image
features with random forest regressor approach was 51.5% for the training and
51.7% for the validation dataset. The gradient boosting regressor with shape
features gave an accuracy of 91.5% and 62.1% on training and validation
datasets respectively. It is better than the BraTS 2020 survival prediction
challenge winners on the training and validation datasets. Our work shows that
handcrafted features exhibit a strong correlation with survival prediction.
The consensus based regressor with gradient boosting and radiomics shape
features is the best combination for survival prediction.

###### keywords:

Brain tumor segmentation (BraTS 2020), glioblastoma, survival prediction

11footnotetext: All authors have contributed equally to this work.

##  1 Introduction

Glioblastoma multiforme (GBM) is the commonest type of primary malignant brain
tumor. In the case of adults, glioblastoma makes up 60% of all brain tumors
[1]. The World Health Organization (WHO) classified GBM as a grade IV type of
cancer due to its invasive and diffusive nature. Patients suffering from GBM
have a poor prognosis, with a median survival rate of about ten months [1].
This is due to its aggressive nature, highly heterogeneous appearance,
location, shape, and unpredictable response to therapy [2].

Magnetic Resonance Imaging (MRI) has been widely utilized to examine tumors
due to its non-hazardousness, high contrast and superior resolution.
Generally, manual segmentation of a tumor in MRI is time consuming and prone
to subjective error. In this regards an automated segmentation method would be
of enormous help to oncologists and clinicians. It can help in early diagnosis
as well as in therapeutic strategy planning. In recent years, deep learning-
based segmentation approaches have outperformed traditional state-of-the-art
methods [3, 4]. Segmentation delineates the brain tumor into Whole Tumor (WT),
Enhancing Tumor (ET), and Tumor Core (TC). Handcrafted features extracted from
these segments are used to classify the survival days of the patients.

There are many segmentation models available. Recently, Jiang et al. [5], in
the BraTS 2019 challenge, proposed a two-stage asymmetry cascaded U-Net [2]
structure. Each model is made up of a larger encoder in order to be able to
extract more complex semantic features and a smaller decoder part for
generating a segmentation map with a size identical to the input. Zhao et al.
[3] proposed multiple methods to generate robust segmentation results. They
grouped it into data processing, model devising, and optimization modules.
Multiple methods are assimilated into each of these modules to enhance
segmentation results. McKinley et al. [4] proposed a Densenet based U-Net
architecture. Convolutions that were dilated were used to bring about an
increase in the receptive field, which retains spatial information. The model
was trained by combining label uncertainty loss, binary cross-entropy and
focal loss. Dice scores on the BraTS-2019 validation dataset were 0.91(WT),
0.83(TC), 0.77(ET), and on the BraTS-2019 test dataset were 0.89(WT),
0.83(TC), 0.81(ET). Therefore, researchers seem to be favouring the U-Net
based architecture for segmentation.

Once the tumor is segmented, features are extracted for overall survival
prediction. Agravat et al. [6] used dense layers U-Net trained on the focal
loss for segmentation. Next, age, statistical features and radiomic features
train the Random Forest Regressor (RFR) for survival prediction and the
obtained accuracy on the test dataset was 0.58. Wang et al. [7] used U-Net and
U-Net ensembles with attention gates trained on soft dice scores and cross-
entropy segmentation. For survival prediction, they proposed the following
prognosis models: i) baseline model where only the age feature was used to
train a linear regressor model. ii) Radiomic model where morphological and
texture features were extracted from segmentation results. iii) Tumor
invasiveness model, where relative invasiveness coefficient (RIC) and age
feature train the support vector regressor model. The tumor invasive model was
found best for survival prediction. The accuracy for survival prediction was
0.59 and 0.56 for BraTS-2019 validation and test dataset respectively. Feng et
al. [8] used an ensemble of U-Net models. The models were trained on patches
having brain pixels. The main advantage of using an ensemble method is that
the network parameter need not be fine-tuned. Further, for OS prediction,
volume and surface area features were extracted for each Region of Interest
(ROIs) and age to train a linear regression model. The training and testing
set accuracy was reported as 0.31 and 0.55 respectively on the BraTS-2019
datasets. Wang et al. [9] utilized a 3D U-Net-based model, and the training
occurred in two phases using patching methods. The first phase included both
brain and background pixels, whereas the second included only brain pixels.
The dice score coefficient loss function was utilized to train the 3D U-Net
model. Further for survival prediction, volume, surface area and age were used
to train the ANN model. The training, validation, and testing accuracy of the
models were 0.515, 0.448, and 0.551 respectively. Islam et al. [10] proposed a
3D U-Net architecture for segmentation, where attention blocks have been
desegregated with the decoder modules. For survival prediction, various
geometric, fractal, and histogram-based features were extracted to train
multiple regressor models, i.e., support vector machine (SVM), multi-layer
perceptron (MLP), random forest regressor (RFR), and eXtreme gradient Boosting
(XGBOOST). The validation accuracies were: 0.329 for SVM, 0.414 for MLP, 0.356
for RFR and 0.429 for XGBOOST.

The proposed paper aims to establish the correlation between handcrafted
features and overall survival prediction. Unlike the existing state-of-the-art
methods used for survival prediction [6],[7], [8], [9], the paper uses four
predictors and two feature sets to establish their correlation with overall
survival prediction of High Grade Glioma (HGG) patients. Shape features and
gradient boosting regressors achieve better survival prediction accuracy than
state-of-the-art methods. It establishes that shape features have a strong
correlation with survival prediction. The organization of the remainder of the
paper is as follows: The Brain Tumor Segmentation (BraTS) dataset is described
in Section 2, survival prediction methods with four predictors and two feature
sets are in Section 3, Section 4 contains results and discussions and finally
the conclusion of the paper is in Section 5.

##  2 BraTS dataset

Due to different standards and differences in the dataset, evaluating brain
tumor segmentation methods objectively and predicting overall survival is a
challenge. Nevertheless, for a comparison of different tumor segmentation and
survival prediction techniques, the BraTS (brain tumor segmentation challenge)
[11, 12, 13] has become a popular platform. Since the year 2018, there are
three tasks that are included in this platform. The first task is the process
of segmenting the brain tumor. The second task is predicting the overall
survival (OS) and the third task is estimating the uncertainty for the
predicted tumor sub-regions. The process of tumor segmentation involves
delineating the tumor into three sub-regions, namely, the whole tumor, the
tumor core, and the enhancing tumor. Specificity and sensitivity metrics as
well as Dice score and Hausdorff Distance are used for evaluating performance.

The overall survival prediction task classifies survival days into the
following categories: long-term survivors ($>$15 months), intermediate-
survivors (between 10 and 15 months), and short-survivors ($<$10 months).
Samples with resection status GTR (gross total resection) are used to rate the
performance of the OS prediction. An accuracy metric is used for performance
evaluation, whereas mean and median square error are used for postanalysis
[14].

The BraTS 2020 training dataset includes 369 volumetric samples of high-grade
glioma (HGG) and low-grade glioma (LGG) cases. It includes metadata of 236
samples such as age, survival days, and resection status for survival days
prediction (Grosstotal Resection (GTR) = 119, Sub-total Resection (STR) = 10,
and NA = 107). The validation dataset includes 125 sample images and metadata
(age, survival days, and resection status) with 29 images having a GTR
resection status. Each subject includes four MRI scans that are preoperative
(T1-weighted, T1-CE, T2-weighted, and FLAIR) and manually annotated ground
truth results. The annotations of ground truth include Necrotic and Non-
Enhancing tumor core NCR/NET (label-1), Edema (label-2), Active Tumor
(label-4), and 0 for everything else. The dataset has been pre-processed,
i.e., all the scans are co-registered to the same anatomical structure, skull
stripped and resampled to an isotropic resolution of $1\times 1\times 1\hskip
5.69046ptmm^{3}$. The width, height, and depth of each sample are 240, 240,
and 155 respectively.

##  3 Survival Prediction Methodology

We use the 3D U-Net model for brain tumor segmentation proposed by Isensee et
al. [15]. This is the highest ranking and simple model in BraTS 2017. Like the
U-Net [2], this model [15] comprises a contracting path to extract more
feature information with increasing network depth. It has an expansion path to
generate a segmentation mask with precise localization information and a skip
connection for better feature reconstruction at every stage of the expansion
path. In our work we have used the bias field correction, normalization,
clipping maximum/ minimum intensity to remove outliers, rescaled to $[0,\,1]$
and setting non-brain pixels to 0. The model was trained on a patch size of
128×128×128, randomly generated from all the input MRI modalities. The
obtained dice score on the BraTS 2020 validation dataset is 0.880(WT),
0.858(TC), 0.759(ET). The segmentation of tumor tissue of a validation sample
is as shown in 1. The figures show a visual comparison of an input flair image
and a predicted image. The segmented parts are then used for survival
prediction with the prognosis methods with 1) Image-based features, 2)
Radiomics based features, and the following four predictors.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 1: Segmentation results of training set: (a) Axial FLAIR slice (b)
Axial Ground truth (c) Axial Segmentation (d) Coronal FLAIR slice (e) Coronal
Ground truth (f) Coronal Segmentation (g) Sagittal FLAIR slice (h) Sagittal
Ground truth (i) Sagittal Segmentation, four color codes are: Brown for
label-1(NCR/NET), white for label-4(Active Tumor), orange for label-2(Edema),
black for label-0(back ground)

###  3.1 Predictors and Parameter Tuning

We have used four predictors and parameter tuning. These are (1) Artificial
Neural Network (ANN) [9, 10], (2) Linear Regressor (LR) [7, 8], (3) Gradient
Boosting Regressor (GBR) [10], and (4) Random Forest Regressor (RFR) [6, 15,
10]. All these predictors were used by the top performing models in all recent
BraTS challenges. These predictors deal with a small dataset and overfitting
problems. The image-based prognosis method uses only seven features making it
less vulnerable to overfitting. We retain default parameters for ANN and LR,
while parameters for GBR and RFR are hyper-tuned using a grid search. We tuned
the number of estimators, depth of the tree, sample split, and learning rate
parameters for the GBR. In the case of the RFR, the number of estimators and
the depth of the tree were hyper tuned. The predictors with radiomics features
were also tuned.

For radiomics features it turns out that an ANN with five hidden layers was
better compared to 2 or 3 hidden layers. Further, we tuned epochs, learning
rate, number of neurons, and an optimizer for ANN. In the LR model, a search
was also performed for the penalty term, the number of iterations, and up-
grading of feature parameters using LASSO and a ridge regressor. We tuned the
number of estimators, maximum depth, and learning rate for the GBR. In the RFR
model, we tuned the number of estimators, maximum depth of the tree, minimum
sample split, minimum samples in a leaf node, and maximum features parameters.
Since the random forest and gradient boosting regressor work on ensemble-based
learning, they are robust, efficient, and less prone to overfitting.

###  3.2 Prognosis using Features

####  3.2.1 Image-based features [8, 9]

Shape features extracted from the segmentation were used in the OS prediction.
These features were volume of the WT, TC, and ET, surface area of the WT, TC,
and ET, age. Since the tumor size was the decisive predicting factor for
various cancer types, we extracted the volume and surface area of the WT, TC,
and ET. The features were extracted from the segmentation maps and input
images without any library dependency. Training with fewer features has the
advantage that it limits the dimensions of feature space. Hence, the model did
not overfit. However, we found saturation in the performance due to high bias
in the model.

####  3.2.2 Radiomics based features [16]

Radiomics based feature extraction is widely used for disease diagnosis,
classification, and survival prediction like lung cancer [17], breast cancer
[18], and Alzheimer’s disease [19]. Along with the size of the tumor,
exploring the correlation of the other features with survival prediction is
crucial to increase the performance of the predictor models. Radiomics
features addresses this problem. It allows extracting various statistical,
shape, intensity, and texture features from radiographic scans. Also,
radiomics allow extracting features from many imaging techniques.

Using the package PyRadiomics [16]: the following 107 features were extracted:

  1. 1.

Shape features: Elongation, flatness, axis lengths, maximum diameter, mesh
volume, sphericity, surface area, and surface volume ratio.

  2. 2.

Gray level features: Gray-level size zone (GLSZ), Gray-level co-occurrence
matrix (GLCM), Gray-level run-length matrix (GLRLM), Gray-level dependence
Matrix (GLDM), and neighbouring gray-tone difference matrix (NGTDM).

  3. 3.

First-order statistical features: Energy, entropy, minimum intensity value,
maximum intensity value, mean, median, Interquartile range, percentiles,
absolute deviation, skewness, variance, kurtosis, and uniformity.

Radiomics features are typically multi-collinear and redundant [20]; hence the
correlation between these features needs to be validated for specific real-
world problems. We performed feature selection through recursive feature
elimination (RFE) [21] to remove weaker features and avoid the curse of
dimensionality. RFE is an example of backward feature elimination. With the
given number of estimators, it selects principal features recursively from the
feature set. It refits the model until the desired number of selected features
is eventually reached. Out of 107 features, we selected 20 best ranking
features.

In summary, the four predictors: ANN, RFR, LR, and GBR, are applied to: i) the
seven image-based features, ii) 107 radiomics features, iii) 20 principal
radiomics features, and iv) only shape radiomics features. Literature [6, 15]
also suggests dominance of shape features so we also used all predictors with
only shape features for survival prediction. We trained the models with all
the resection status (i.e., GTR, STR, and NA) given with the dataset to
increase the database size and reduce overfitting.

##  4 Results and Discussions

Image-based feature prediction is derived from the BraTS 2019 dataset, and the
BraTS 2020 dataset was used for radiomics based feature extraction. The
results are shown in Tables 1 to 4. We have not participated in the BraTS 2020
challenge and do not have access to the test dataset. Therefore, results are
derived on the training and validation datasets.

###  4.1 Image-based feature prediction

We observe that the ensemble-based models, i.e., GBR and RFR, show a better
performance on the training and validation dataset. Their consistency in the
training and validation accuracy suggests that the model does not overfit.

Table 1: OS Performance comparison using image-based feature on training and validation BraTS-2019 dataset. MSE, medianSE, stdSE, and SpearmanR denote the mean square error, median square error, standard deviation squared error, and Spearman’s ranking coefficient. Dataset | Regressor | Accuracy | MSE | medianSE | stdSE | SpearmanR
---|---|---|---|---|---|---
Training | ANN | 0.51 | 86148.10 | 21316 | 181346 | 0.48
LR | 0.49 | 87724.00 | 20736 | 183685 | 0.47
GBR | 0.52 | 63234.40 | 16900 | 126534 | 0.61
RFR | 0.52 | 63234.40 | 16900 | 126534 | 0.61
Validation | ANN | 0.45 | 098312.70 | 39204 | 141392 | 0.24
LR | 0.52 | 100509.00 | 38809 | 141263 | 0.29
GBR | 0.52 | 102999.00 | 36481 | 152694 | 0.27
RFR | 0.52 | 102999.00 | 36481 | 152694 | 0.27

###  4.2 Radiomics feature-based prediction

As mentioned, we extracted 107 radiomic features from the segmentation results
of the BraTS 2020 images and fed them as input to four regressor models; ANN,
LR, GBR, and RFR. It was observed that RFR gave the best results, and they are
shown in Table 2. The other regressors performed poorly compared to RFR, and
even the fine-tuning of the parameters did not improve the performance. The
possible reasons are the redundant nature of radiomics [20], over complexity
due to too many features and fewer training samples. Radiomics features are
shallow and low-order image features, and unable to fully describe distinct
image characteristics [22]. Also, when the number of observations is less for
large extracted features, survival prediction is an ill-posed problem [20].

Table 2: OS performance evaluation using 107 radiomics features and Random Forest Regressor.  Dataset |  Accuracy |  MSE |  medianSE |  stdSE |  SpearmanR
---|---|---|---|---|---
Training |  0.479 |  079176.96 |  20702.21 |  169474.53 |  0.684
Validation |  0.379 |  115424.30 |  28779.30 |  214028.11 |  0.138

It can be observed from Table 2 that the large feature set is unable to yield
state-of-the-art accuracy results. Therefore, we reduced the feature set by
applying recursive feature elimination to find the 20 most dominant features.
Dominant features obtained using RFE are: age, amount of edema, elongation,
maximum 2D diameter slice, sphericity, surface-volume ratio, minimum and
maximum intensity, interquartile range, skewness, kurtosis, root mean absolute
deviation, cluster prominence, cluster shade, inverse variance, coarseness,
and dependence variance. We then applied four regressors on the dominant
feature set, and performance has been noted in Table 3.

Table 3: OS performance comparison on 20 principal radiomics features. Dataset |  Regressor Models |  Accuracy |  MSE |  medianSE |  stdSE |  SpearmanR
---|---|---|---|---|---|---
Training |  ANN |  0.393 |  8.90E+12 |  2.46E+12 |  3.36E+13 |  0.125
LR |  0.462 |  96853.55 |  33279.52 |  190733.00 |  0.417
GBR |  0.923 |  17213.25 |  00000.00 |  074717.13 |  0.938
RFR |  0.744 |  31829.75 |  06077.32 |  075572.44 |  0.810
Validation |  ANN |  0.448 |  2.20E+20 |  3.46E+12 |  8.03E+20 |  0.290
LR |  0.483 |  2.73E+08 |  056167.55 |  9.86E+08 |  0.456
GBR |  0.414 |  255096.40 |  101995.06 |  420861.25 |  0.025
RFR |  0.448 |  098369.46 |  035521.48 |  126218.18 |  0.126

We observe that the linear regressor with regularisation outperforms all other
regression models with the highest accuracy on the validation dataset. LR also
provides similar accuracy for the training and validation datasets. The
Spearman-R is also highest for LR. In contrast, RFR achieves the lowest mean
square error (MSE) on the validation dataset.

####  4.2.1 Radiomic shape features based prediction

Reviewing the correlation between radiomics features and survival prediction,
we found that radiomic shape features play a crucial role in survival
prediction [6, 15]. Shape features show significant statistical differences
across ROIs [23]. Hence, shape features can capture tumor features related to
genetic anomalies and profoundly impact survival prediction. We formulate the
hypothesis that _shape features profoundly impact survival prediction_. In
order to validate the hypothesis, we trained predictor models with the
following shape features: the amount of necrotic, edema, enhancing tumor, the
extent of the tumor, coordinates of tumor, elongation, flatness, axis lengths,
2D diameter row, 2D diameter column, 2 D diameter slice, maximum 3D diameter,
mesh volume, sphericity, surface area, surface volume ratio, centroid of
necrosis and age information. The performance of each predictor model has been
noted in Table 4.

Table 4: OS performance comparison on BraTS-2020 dataset using radiomics shape features set. Dataset |  Predictor Models |  Accuracy |  MSE |  medianSE |  stdSE |  SpearmanR
---|---|---|---|---|---|---
Training |  ANN |  0.400 |  4.41E+11 |  7.15E+10 |  7.97E+11 |  0.149
LR |  0.470 |  89890.41 |  35160.09 |  162137.20 |  0.461
GBR |  0.915 |  31068.75 |  00000.00 |  150724.63 |  0.849
RFR |  0.615 |  62930.78 |  18562.88 |  130788.18 |  0.759
Validation |  ANN |  0.448 |  4.73E+11 |  2.14E+11 |  5.97E+11 |  0.149
LR |  0.414 |  087228.24 |  47820.00 |  111960.30 |  0.215
GBR |  0.621 |  141065.30 |  23528.48 |  236728.70 |  0.338
RFR |  0.448 |  109746.60 |  34689.29 |  200725.98 |  0.116

We observe that GBR and RFR have better performance. Specifically, the
gradient boosting regressor outperforms all other regression models. In
contrast, LR with regularization achieves the lowest mean square error (MSE)
on the validation dataset.

###  4.3 Discussions

It has been observed that classical machine learning techniques performed
better than the deep learning neural network-based models for survival
prediction. Radiomics based approaches are well suited for survival
prediction. Traditional regression algorithms have better interpretability
than deep learning-based algorithms, they have fewer learnable parameters than
CNN, and perform better with smaller sample data. A large sample dataset for
training is crucial for direct regression from image modalities using CNN.

The predictors trained on the 107 radiomics features underperformed. The
predictors modelled on the 20 principal features improved the performance.
Further, to alleviate performance, we experimented and trained predictors on
shape features and found a strong correlation with survival prediction. Shape
features trained on the consensus model obtained state-of-the-art survival
prediction accuracy. It was observed that the gradient boosting regressor
model performed better than other classical algorithms because of: additive
model, and with each tree built, the model becomes more expressive based on
the ensemble learning model. The proposed GBR model is compared with the
survival prediction challenge winners of BraTS 2020 and prediction accuracy
for the state-of-the-art methods was obtained from the unranked leader
board***https://www.cbica.upenn.edu/BraTS20/lboardValidation.html. TA
performance comparison of the GBR model with top-ranking models has been noted
in Table 5. It can be observed that shape-based features with the gradient
boosting regressor outperform the best-ranking methods over the validation
dataset.

Table 5: OS performance comparison with top-ranking models on the BraTS-2020 validation dataset. Team name |  Accuracy |  MSE |  medianSE |  stdSE |  SpearmanR
---|---|---|---|---|---
SCAN |  0.414 |  098704.65 |  36100.00 |  152175.57 |  0.253
Redneucon |  0.517 |  122515.76 |  70305.26 |  157673.99 |  0.134
VLB |  0.379 |  093859.54 |  67348.26 |  102092.41 |  0.280
COMSATS-MIDL |  0.483 |  105079.42 |  37004.93 |  146375.99 |  0.134
Proposed |  0.621 |  141065.30 |  23528.40 |  236728.70 |  0.338

##  5 Conclusion

Predicting oncological outcomes is always very tricky due to multiple
challenges from clinical and engineering perspectives. In this work, we have
evaluated two feature sets over four predictors. We proposed the image-based
and the radiomic based prognosis approaches for survival prediction. The
image-based prognosis models performed well, but the performance saturates
beyond a certain point because of fewer features, and models could not learn
complexity. Similar observations are also made for the 107 radiomics features
/ 20 principal features and the regressor combination. All above the
combinations exhibited correlation with survival prediction. However, we
recommend that shape based features with the gradient boosting regressor is
the best combination for survival prediction. Comparing models, it was found
that ensemble-based learning models became more useful for survival prediction
because of their robustness. Whereas ANN converges speedily compared to
classical models but due to lack of ample training samples, it overfits
easily. With the availability of a large dataset and more clinical non-imaging
information such as gender and treatment, survival prediction can be robust.
It can further be applied to clinical practice.

## References

  * [1] Taylor, O.G., Brzozowski, J.S., Skelding, K.A.: Glioblastoma multiforme: An overview of emerging therapeutic targets. Frontiers in oncology 9, 963 (2019)
  * [2] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [3] Zhao, Y.X., Zhang, Y.M., Liu, C.L.: Bag of tricks for 3d mri brain tumor segmentation. In: International MICCAI Brainlesion Workshop. pp. 210–220. Springer (2019)
  * [4] McKinley, R., Rebsamen, M., Meier, R., Wiest, R.: Triplanar ensemble of 3d-to-2d cnns with label-uncertainty for brain tumor segmentation. In: International MICCAI Brainlesion Workshop. pp. 379–387. Springer (2019)
  * [5] Jiang, Z., Ding, C., Liu, M., Tao, D.: Two-stage cascaded u-net: 1st place solution to brats challenge 2019 segmentation task. In: International MICCAI Brainlesion Workshop. pp. 231–241. Springer (2019)
  * [6] Agravat, R.R., Raval, M.S.: Brain tumor segmentation and survival prediction. In: International MICCAI Brainlesion Workshop. pp. 338–348. Springer (2019)
  * [7] Wang, S., Dai, C., Mo, Y., Angelini, E., Guo, Y., Bai, W.: Automatic brain tumour segmentation and biophysics-guided survival prediction. In: International MICCAI Brainlesion Workshop. pp. 61–72. Springer (2019)
  * [8] Feng, X., Dou, Q., Tustison, N., Meyer, C.: Brain tumor segmentation with uncertainty estimation and overall survival prediction. In: International MICCAI Brainlesion Workshop. pp. 304–314. Springer (2019)
  * [9] Wang, F., Jiang, R., Zheng, L., Meng, C., Biswal, B.: 3d u-net based brain tumor segmentation and survival days prediction. In: International MICCAI Brainlesion Workshop. pp. 131–141. Springer (2019)
  * [10] Islam, M., Vibashan, V., Jose, V.J.M., Wijethilake, N., Utkarsh, U., Ren, H.: Brain tumor segmentation and survival prediction using 3d attention unet. In: International MICCAI Brainlesion Workshop. pp. 262–272. Springer (2019)
  * [11] Bakas, S., Akbari, H., Sotiras, A., Bilello, M., Rozycki, M., Kirby, J.S., Freymann, J.B., Farahani, K., Davatzikos, C.: Advancing the cancer genome atlas glioma mri collections with expert segmentation labels and radiomic features. Scientific data 4, 170117 (2017)
  * [12] Bakas, S., Reyes, M., Jakab, A., Bauer, S., Rempfler, M., Crimi, A., Shinohara, R.T., Berger, C., Ha, S.M., Rozycki, M., et al.: Identifying the best machine learning algorithms for brain tumor segmentation, progression assessment, and overall survival prediction in the brats challenge. arXiv preprint arXiv:1811.02629 (2018)
  * [13] Menze, B.H., Jakab, A., Bauer, S., Kalpathy-Cramer, J., Farahani, K., Kirby, J., Burren, Y., Porz, N., Slotboom, J., Wiest, R., et al.: The multimodal brain tumor image segmentation benchmark (brats). IEEE transactions on medical imaging 34(10), 1993–2024 (2014)
  * [14] Rajput, S., Raval, M.S.: A review on end-to-end methods for brain tumor segmentation and overall survival prediction. arXiv preprint arXiv:2006.01632 (2020)
  * [15] Isensee, F., Kickingereder, P., Wick, W., Bendszus, M., Maier-Hein, K.H.: Brain tumor segmentation and radiomics survival prediction: Contribution to the brats 2017 challenge. In: International MICCAI Brainlesion Workshop. pp. 287–297. Springer (2017)
  * [16] Van Griethuysen, J.J., Fedorov, A., Parmar, C., Hosny, A., Aucoin, N., Narayan, V., Beets-Tan, R.G., Fillion-Robin, J.C., Pieper, S., Aerts, H.J.: Computational radiomics system to decode the radiographic phenotype. Cancer research 77(21), e104–e107 (2017)
  * [17] He, B., Zhao, W., Pi, J.Y., Han, D., Jiang, Y.M., Zhang, Z.G.: A biomarker basing on radiomics for the prediction of overall survival in non–small cell lung cancer patients. Respiratory research 19(1), 1–8 (2018)
  * [18] Liu, C., Ding, J., Spuhler, K., Gao, Y., Serrano Sosa, M., Moriarty, M., Hussain, S., He, X., Liang, C., Huang, C.: Preoperative prediction of sentinel lymph node metastasis in breast cancer by radiomic signatures from dynamic contrast-enhanced mri. Journal of Magnetic Resonance Imaging 49(1), 131–140 (2019)
  * [19] Li, Y., Jiang, J., Lu, J., Jiang, J., Zhang, H., Zuo, C.: Radiomics: a novel feature extraction method for brain neuron degeneration disease using 18f-fdg pet imaging and its implementation for alzheimer’s disease and mild cognitive impairment. Therapeutic Advances in Neurological Disorders 12, 1756286419838682 (2019)
  * [20] Weninger, L., Haarburger, C., Merhof, D.: Robustness of radiomics for survival prediction of brain tumor patients depending on resection status. Frontiers in computational neuroscience 13, 73 (2019)
  * [21] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al.: Scikit-learn: Machine learning in python. the Journal of machine Learning research 12, 2825–2830 (2011)
  * [22] Lao, J., Chen, Y., Li, Z.C., Li, Q., Zhang, J., Liu, J., Zhai, G.: A deep learning-based radiomics model for prediction of survival in glioblastoma multiforme. Scientific reports 7(1), 1–8 (2017)
  * [23] Chaddad, A., Desrosiers, C., Hassan, L., Tanougast, C.: A quantitative study of shape descriptors from glioblastoma multiforme phenotypes for predicting survival outcome. The British journal of radiology 89(1068), 20160575 (2016)


# Interplay of Andreev reflection and Coulomb blockade in hybrid
superconducting single electron transistors

Laura Sobral Rey  Physics Department, University of Konstanz, 78457 Konstanz,
Germany  David Christian Ohnmacht  Physics Department, University of Konstanz,
78457 Konstanz, Germany  Clemens B. Winkelmann  Univ. Grenoble Alpes, CNRS,
Grenoble INP, Institut Néel, 38000 Grenoble, France  Jens Siewert  University
of the Basque Country UPV/EHU and EHU Quantum Center, 48080 Bilbao, Spain
Ikerbasque, Basque Foundation for Science, 48009 Bilbao, Spain  Wolfgang
Belzig  Physics Department, University of Konstanz, 78457 Konstanz, Germany
Elke Scheer<EMAIL_ADDRESS>Physics Department, University of
Konstanz, 78457 Konstanz, Germany

(August 28, 2024)

###### Abstract

We study the interplay between Coulomb blockade and superconductivity in a
tunable superconductor-superconductor-normal metal single-electron transistor.
The device is realized by connecting the superconducting island via an oxide
barrier to the normal metal lead and with a break junction to the
superconducting lead. The latter enables Cooper pair transport and (multiple)
Andreev reflection. We show that those processes are relevant also far above
the superconducting gap and that signatures of Coulomb blockade may reoccur at
high bias while they are absent for small bias in the strong-coupling regime.
Our experimental findings agree with simulations using a master equation
approach in combination with the full counting statistics of multiple Andreev
reflection.

Coulomb blockade (CB) is an archetypical manifestation of charge quantization
(CQ), which occurs in the electronic transport across a small metallic island
[1]. CB can be suppressed by both classical and quantum fluctuations of the
charge [2]. Classical charge fluctuations originate from thermal activation
over an energy barrier, provided by the charging energy $E_{\rm C}=e^{2}/2C$,
with $C$ the total capacitance and $e$ the elementary charge. The quantum
fluctuations can be described using the Landauer-Büttiker picture using
individual conductance channels $i=1\dots N$, with transmissions
$0\leq\tau_{i}\leq 1$, connecting the island to the leads [2]. In the case of
normal metallic leads, the magnitude of the conductance oscillations in a
Coulomb-blockaded island in the single electron transistor (SET) geometry [3,
4] was predicted to scale like
$\prod_{i}\sqrt{1-\tau_{i}}\exp(-k_{\text{B}}T/E_{\text{C}})$ [5], which was
shown to hold experimentally with great accuracy [6]. In the presence of
superconducting (S) contacts this picture can be expected to change radically
due to the different nature of the charge carriers.

For mesoscopic transport processes involving S leads, the energy-dependent
quasiparticle (QP) spectrum has a gap $\Delta$, which strongly inhibits
tunneling at low energies and small $\tau$ [7]. At larger $\tau$,
multiparticle superconductive transport (MST) [7, 8, 9] comes into play. These
processes can be either coherent, in the form of Josephson transport of Cooper
pairs (CPs), in the absence of a voltage drop, or they may be dissipative, as,
e.g., multiple Andreev reflection (MAR) of order $m$, which sets in above a
bias voltage threshold $eV=2\Delta/m$. Importantly, while MAR is itself
disspative and therefore incoherent, the transfer of $m$ charges involved in
such a MAR is a coherent tunneling process. The theoretical description and
to-date experimental realizations of MAR were based on superconducting
contacts described by a well-defined quantum phase, i.e., large particle
number fluctuations.

From the above discussion, it is clear that MST and CQ are antagonistic
processes. On the one hand, CQ effects should deeply modify the energy
thresholds for a $m$-particle tunneling process On the other hand, MST could
overthrow the CQ conditions known from the case of normal leads.

Figure 1: Scanning electron microscope picture of a single electron transistor
(SET) with a break junction (BJ) as one of the tunnel junctions in series with
a superconductor-insulator-normal metal tunnel junction (TB), and the area in
the middle (island) being capacitively coupled to a gate electrode.

In this Letter, we address the delicate interplay of CQ and MST. To this end
we build a superconductor-superconductor-normal metal (SSN)-SET, hosting a
superconducting break junction (BJ), with a small set of mechanically tunable
transmission channels. The BJ is connected in series with a superconductor-
normal metal (SN) tunnel barrier (TB), providing thus the conditions for CQ in
the S island formed between the BJ and the TB, see Fig. 1. The device is
completed by a capacitively coupled gate electrode [3, 4]. In the limit of a
small total transmission
$\mathcal{T}_{\text{BJ}}=\Sigma_{i}\tau_{i,\text{BJ}}\ll 1$ of the BJ, we
study how the extra energy cost arising from CB affects the conductance
thresholds for the different MAR processes. Furthermore, in the limit of a
larger $\mathcal{T}_{\text{BJ}}\sim 1$, we observe that CQ enters in
competition with the Josephson effect across the BJ, leading to the complete
suppression of CB at low bias, whereas it is entirely restored beyond the
Josephson critical current. Our data are well described by a master equation
approach including CB as well as MST.

The transport in the SET is characterized by current cycles composed of charge
transfer processes through the individual junctions comprising different types
and different number of charged particles [10]. The $e$-cycle is a process
that charges the island with one QP through one junction, and discharges it
through the other junction with another QP, leaving the island in its initial
charge state. AR and MAR cycles charge the island with $m$ electron charges
through one junction, and subsequent single-QP or (M)AR processes discharge
the island through the other junction. The Josephson quasiparticle (JQP) cycle
[11, 12, 13, 14] is a coherent process in which a CP tunnels through one
junction, and simultaneously one QP tunnels off through the other junction.
The initial charge state is reached by a subsequent single-QP process. The JQP
cycle is resonant, i.e., it appears as a peak in the current-bias voltage
($I-V_{\text{b}}$) characteristics.

Studies on all-superconducting (SSS)-SETs in different realizations were
conducted in the weak-coupling regime in which both junctions were realized as
tunnel contacts and $e$, JQP and processes comprising 3 QPs were found [10].
Another realization using a tunable junction showed a multitude of current
cycles which were difficult to identify [15]. Hence the SSN-SET studied here
is chosen to enable MST predominantly through just one of the junctions (BJ),
thereby reducing the number of processes. While some theoretical studies on
the SSN-SET exist [16, 17], to the best of our knowledge, it has never been
investigated experimentally before.

For forming the tunable junction we use the thin-film BJ technique [18]. The
BJ consists of an Al bridge suspended above a flexible substrate: by bending
the substrate it is possible to adjust contacts with arbitrary transmission
$\mathcal{T}_{\text{BJ}}$ from the tunnel regime $\mathcal{T}_{\text{BJ}}\ll
1$ to atomic-size contacts with $\mathcal{T}_{\text{BJ}}\gtrsim 1$ [19, 20].
The TB with fixed resistance $R_{\text{TB}}$ is formed by an AlxOy layer
between the S (Al) and the N (Cu) material and hosts many channels in parallel
with very small transmission $\tau_{i}\ll 1$ [2]. The electrostatic energy
necessary to charge the island is $\mathcal{E}(n_{\rm
g},n)=E_{\text{C}}(n-n_{\text{g}})^{2}$ where $n$ is the excess charge (in
units of $e$) of the island, $n_{\text{g}}$ is the equivalent charge induced
by the gate voltage, and $C=C_{\text{BJ}}+C_{\text{TB}}+C_{\text{g}}$ the
total capacitance of the island. Furthermore we define the capacitance
divisors $\kappa_{i}=(C_{i}+C_{\text{g}}/2)/C$. In this Letter, we discuss
data taken on one sample.The sample preparation and parameters as well as data
on three other samples with similar properties are given in the Supplemental
Material (SM) [21].

Figure 2: Evolution of the transport processes in the SSN-SET for different
$\mathcal{T}_{\text{BJ}}$: (a) $dI/dV_{\text{b}}-eV_{\text{b}}$ cuts for
varying $\mathcal{T}_{\text{BJ}}$ as indicated in the legend.
$n_{\text{g}}-eV_{\text{b}}-dI/dV_{\text{b}}$ maps for
$\mathcal{T}_{\text{BJ}}=0.01$ (b) and $\mathcal{T}_{\text{BJ}}=3$ (c).

We start the discussion by comparing different coupling strengths for the BJ,
from weak to strong. Fig. 2(a) shows $dI/dV_{\text{b}}-eV_{\text{b}}$ plots
taken at fixed $n_{\text{g}}$ for different $\mathcal{T}_{\text{BJ}}$ (from
top to bottom) of 0.01, 0.07, 0.4, 0.72 and 3. From these individual
conductance curves taken at different $n_{\text{g}}$ we construct conductance
vs. gate voltage and bias voltage maps in reduced units as
$n_{\text{g}}-V_{\text{b}}-dI/dV_{\text{b}}$ color plots. Two extreme
examples, for $\mathcal{T}_{\text{BJ}}=0.01$ and 3, are given in Fig. 2(b,c),
respectively. For the other transmissions, we plot the
$n_{\text{g}}-eV_{\text{b}}-I$ in Fig. 3(a,b) ($\mathcal{T}_{\text{BJ}}=0.07$
and 0.4, respectively) and Fig. 4 for $\mathcal{T}_{\text{BJ}}=0.72$. Further
examples are shown in the SM [21]. The maps show that different current
contributions are modulated by the gate voltage with $e$-periodicity with
asymmetric slopes indicating different capacitances of the junctions. The
positive/negative slope mark the onset of charge transport across the TB/BJ,
respectively. The visibility of the transport across the BJ is much higher
than the one across the TB as a consequence of the different nature of the
junctions. The vertical red lines mark the positions at which the data in Fig.
2(a) have been taken. They correspond to positions in the middle of the
descending slope of a CB diamond 111We show here data for positive
$V_{\text{b}}$, but note that all maps were symmetric with respect to the bias
polarity when choosing the adequate $n_{\rm g}$. For intermediate
$\mathcal{T}_{\text{BJ}}\sim 0.5$ the peak slightly above $\Delta$ is
hysteretic with respect to the bias sweep direction.. For more details see the
SM [21].

For the lowest $\mathcal{T}_{\text{BJ}}=0.01$ in Fig. 2(a), we observe only
one clearly visible peak at $3\Delta$, marking the $e$-process. For
$\mathcal{T}_{\text{BJ}}=0.07$, the peak corresponding to the $e$-process is
bigger in amplitude, and two additional peaks appear between $\Delta$ and
$3\Delta$. The first, that – depending on $n_{\rm g}$ – is located between
$\Delta+E_{\text{C}}$ and $\Delta+3E_{\text{C}}$ ($\sim$1.4$\Delta$ and
2.2$\Delta$), corresponds to a JQP cycle [15, 10, 23, 13], while we assign the
second one (with a lower threshold of $2\Delta+E_{\text{C}}$) to a cycle
involving AR across the BJ. Finally, the peak at around $4.5\Delta$ signals
again the next step in the Coulomb staircase. For
$\mathcal{T}_{\text{BJ}}=0.4$ the amplitude of the JQP and the AR cycle
increase further. At the position of the $e$-process a double peak occurs that
we identify as an AR process followed by a MAR process with $m=3$, indicating
that both AR and MAR are relevant processes in this coupling regime. For
$\mathcal{T}_{\text{BJ}}=0.72$ the amplitude of the MAR and JQP cycles has
further increased such that the JQP of the adjacent diamond, i.e., signaling
the next charge state, gives rise to the highest peak at around 2.4$\Delta$.
There is no particular peak starting at 3$\Delta$, meaning that cycles
including MST are more relevant than $e$ transport. Besides, there is a small
peak slightly above $\Delta$ that can be explained by considering cotunneling
of several QPs, see below.For large $\mathcal{T}_{\text{BJ}}$ large Josephson
currents can flow without a voltage drop. In this situationthe voltage drops
completely across the SN-TB, and the island’s chemical potential remains
pinned to that of the S contact. This is visible in the lowest trace in Fig.
2(a) as well as in Fig. 2(c), where above $\Delta/e$ no gate dependence is
observed, indicating a suppression of CQ. Upon further increasing $|eV_{\rm
b}|$ to about $2.2\Delta$, gate-dependent conductance features with the same
slope as observed for small $\mathcal{T}_{\text{BJ}}$ reappear. At this
threshold, the current amounts to 7.0 nA, an order of magnitude below the
estimation $I_{\text{c}}=\frac{\pi}{2}G_{0}\tau\Delta\approx 69$ nA [24] of
the maximally possible supercurrent, in agreement with reports from earlier
experimental [25, 26] and theoretical works [27] in SSS-SETs. As the
dissipation-free state of the SS junction is lost, the island’s chemical
potential is no longer pinned via the Josephson effect to that of the S lead.
Consequently, the SS junction switches from a vanishing to a high
impedance.This restores CB effects despite of the large
$\mathcal{T}_{\text{BJ}}\approx 3$.

Figure 3: Theoretical and experimental $n_{\text{g}}-eV_{\text{b}}-I$ maps for
a contact with $\mathcal{T}_{\text{BJ}}=$ 0.07 (a) and
$\mathcal{T}_{\text{BJ}}$ = 0.4 (b). Red lines in (a) and (b): positions of
cuts of Fig. 2. (c,d): correspondent $I-eV_{\text{b}}$ cuts at constant
$n_{\text{g}}$as marked in (a) and (b) with solid black lines.Black/gray:
experiment (solid/dashed line: increasing/decreasing $|V_{\text{b}}|$)
/simulation. Scatter plots: $P_{n}\Gamma_{n\rightarrow n+n^{\prime}}$
contributions (with amplitude $>0.05e\Delta/h$): $e$ ($e$: $n\rightarrow n+1$,
green), AR (A: $n\rightarrow n+2$, pink), $m=3$ MAR (M3: $n\rightarrow n+3$,
purple), $m=4$ MAR (M4: $n\rightarrow n+4$, yellow). . Insert in (d):
experimental (black) and simulated (gray) $dI/dV_{\text{b}}$.

To verify this assignment of the various contributions, we calculate the
current through the device with the master equation approach (Eq. 2) [1]. We
take into account the JQP rate $\Gamma_{\text{JQP}}(\delta E)$ from Averin and
Aleshkin [28] and rates derived from the full counting statistics for (M)AR
processes $\Gamma_{\text{MAR}}(\delta E)$ [29] in the BJ. For the TB, we
solely use single-QP transport using the classical rate expression for a SN
junction [7] 222For describing the cotunneling processes at intermediate
coupling we will also take a small AR contribution through the TB junction
into account.. The electrostatic energy difference can be expressed as [1]:

$\delta
E=\mathcal{E}(n_{\text{g}},n)-\mathcal{E}(n_{\text{g}},n+n^{\prime})+n^{\prime}\kappa_{i}eV_{\text{b}}$
(1)

where $n^{\prime}$ is the number of transferred charges. We assume that there
is no accumulation of electrons in the island or relevant environment and
neglect effects of the environment [31], so $dP/dt=0$ and the current through
the two junctions is equal. Considering the normalization condition
$\sum_{n}P_{n}=1$ we solve the stationary master equation Eq. (2) and obtain:

$\sum_{n^{\prime}\neq 0}P_{n+n^{\prime}}\Gamma_{n+n^{\prime}\rightarrow n}\ -\
P_{n}\sum_{n^{\prime}\neq 0}\Gamma_{n\rightarrow n+n^{\prime}}\ =\ 0\ \ .$ (2)

From these, we calculate the total current through the SSN-SET by considering
the charge transferred, e.g., in the BJ,

$\begin{split}I\ =\ (-e)\sum_{n;n^{\prime}\neq 0}\tilde{n}\ P_{n}\
\Gamma_{\mathrm{BJ},n\rightarrow n+n^{\prime}}\ \ .\end{split}$ (3)

Here, $n^{\prime}$ can also be negative, and $\tilde{n}$ is the charge
transferred in the BJ (which may be different from $n^{\prime}$).

Figs. 3(a) and (c) are the $n_{\text{g}}-eV_{\text{b}}-I$ map and a
$I-eV_{\text{b}}$ cut at $n_{\text{g}}$ = 0.75 as indicated in the map, for a
configuration with a $\mathcal{T}_{\text{BJ}}=0.07$ for the BJ. In Fig. 3(c),
the black lines correspond to experimental curves (solid/dotted for
increasing/decreasing voltage) while the gray lines correspond to the
simulation. The colored symbols depict the individual contributions
$P_{n}\Gamma_{n\rightarrow n+n^{\prime}}$ that sum up to the total currentThe
resonant peak measured in the experiment follows the JQP contribution from the
simulation, and the small feature starting at 2.4$\Delta$ is due to a cycle
with AR [10, 32]. Fig. 3(b,d) show the comparison of experiment and theory for
a higher transmission $\mathcal{T}_{\text{BJ}}=0.4$. As Fig. 3(b) shows, the
AR, MAR and JQP cycles become more pronounced and the $e$-cycle starting at
$3\Delta$ is almost invisible. The large blue area in the experimental map
above the JQP line (at $\sim 1.5\Delta$) is a result of the resonant character
of the JQP and is absent for the downsweep (see SM). The decomposition into
individual cycles in Fig. 3(d) confirms that AR and MAR with $m=3$ (pink and
purple symbols) are bigger in amplitude than the $e$ transport (green
symbols). In single SS junctions $e$ and AR persist to very high bias, while
MAR fades out above $2\Delta$ [29]. In contrast, here, also MAR contributes to
the current far outside the gap region. This finding shows that MAR processes
become more robust by the interplay with CB. The $dI/dV_{\text{b}}$ in the
inset shows that the multipeak characteristic of each charge state found in
the simulation, indicative for the presence of the MAR processes, is also
observed in the experimental curve 333We note that Al single-atom contacts
with $\mathcal{T}\gtrsim 0.8$ are expected to accommodate more than one
channel [19]. For this contact with $\mathcal{T}=0.4$, it was sufficient to
consider one single channel with $\tau=0.4$. Considering more than one channel
would result in a smaller current, and hence, would not improve the agreement
between experiment and simulation giving a worse agreement..

Figure 4: Experimental and simulated maps for a contact with
$\mathcal{T}_{\text{BJ}}$ = 0.72. The red line marks the position at which the
line cut shown in Fig. 2(a) has been taken.

For even higher $\mathcal{T}_{\text{BJ}}$ (Fig. 4) the amplitudes of the
processes setting in below $3\Delta$ increase, and besides a new line appears
at $|eV_{\text{b}}|=\Delta$. In the $n_{\text{g}}-eV_{\text{b}}$ map the slope
of this cycle is around 3 times smaller than the one of the other cycles.

Here, the lowest-order cycle with the observed slope consists of a coherent
tunneling process that simultaneously transports three charges through the BJ
and two through the TB, hence implies AR in the TB. This process has the
following electrostatic energy balance $\delta
E=\mathcal{E}(n_{\text{g}},n)-\mathcal{E}(n_{\text{g}},n+1)+(3\kappa_{\text{BJ}}+2\kappa_{\text{TB}})eV_{\text{b}}$.
This equation implies that the island charge changes by $e$ and shows that the
slope in the $n_{\text{g}}-eV_{\text{b}}$ map is a function of the capacitance
ratio and the number of charge carriers involved [34, 35]. For the sample
studied here it turns out to be 1/2.88, hence close to 1/3. The expression for
such a rate has not been worked out theoretically. Therefore, we tentatively
assume the rate as a step function starting at $\delta E=2\Delta$ with an
amplitude as free parameter. In Fig. 4 we show the comparison between
experiment and theory for a contact with $\mathcal{T}_{\text{BJ}}$ = 0.72 that
shows this feature starting at $|eV_{\text{b}}|=\Delta$ and with a slope in
agreement with the experimental result. The contributions at higher bias,
namely the JQP and MAR cycles, appear in the simulation at around the same
position as in the experimental map and they follow the same slopes. The
amplitude of both processes overcomes the $e$ transport for high voltages and
extends to very high $|V_{\text{b}}|$. The apparent double-line pattern above
$2.5\Delta$ is the result of JQP and MAR cycles from different initial charge
states of the island. We note that we also included the three-QP cotunneling
process (AR in the BJ and $e$ in the TB) in the simulation. This contribution
sets in at $eV_{\text{b}}\sim 1.5\Delta$, has a slope of $\sim 1/2$ and its
onset is therefore hidden by the JQP contribution. However, it contributes to
the overall current level, improving the agreement between theory and
experiment. Simulation results without these contributions are shown in the SM
[21].

In conclusion, we presented transport properties of an SSN-SET with a tunable
junction allowing to vary the relative strength of $e$, JQP, AR and MAR
processes. With this approach we show the mutual influence between MST
processes and CB. While in single junctions MAR processes are not relevant
outside the gap region and only $e$ and AR processes make up the excess
current, the combination with CB enhances the contributions of MAR at high
bias. We show that the description based on charge states of the island allows
to quantitatively reproduce the experimental data up to relatively strong
coupling $\mathcal{T}_{\text{BJ}}\simeq 1$ when taking all four kinds of
cycles into account. In particular MAR cycles show up in three different ways
in the SET: MAR of order $m$ charges the island with $me$ and its amplitude
can be higher than the one for $e$ transport. AR and MAR are the dominating
transport process up to very high voltages. For higher
$\mathcal{T}_{\text{BJ}}$, cotunneling events become relevant and it is
possible to identify, qualitatively, a coherent multiparticle cotunneling
process across both junctions inside the gap region $|eV_{\text{b}}|\leq
3\Delta$, which implies MAR with $m=3$ in the BJ and AR in the TB. The third
and most striking signature of MAR is the reappearance of charging effects at
high bias in the strong coupling regime ($\mathcal{T}_{\text{BJ}}>1$) when at
least one channel with a transmission $\tau=1$ exists in the BJ) while at low
bias all signature of CQ has already disappeared. Our findings are important
for designing and understanding nanoscale superconducting quantum devices,
such as superconducting qubits [36], with arbitrary intermediate transmissions
to use the subtle interplay between the charge and phase degrees of freedom as
a resource for optimizing their performance.

Acknowledgements: We thank J. C. Cuevas for fruitful discussion and S.
Sprenger and T. Lorenz-Sprenger for experimental support. We gratefully
acknowledge funding by the European Commission through the Marie Skłodowska-
Curie grant agreement no. 766025 (QuESTech) and the Deutsche
Forschungsgemeinschaft (DFG; German Research Foundation) by SFB 1432 (Project
No. 425217212).

## References

  * Grabert and Devoret [1992] H. Grabert and M. H. Devoret, eds.,  _Single Charge Tunneling_ (Springer, Boston, MA, 1992).
  * Nazarov and Blanter [2009] Y. Nazarov and Y. Blanter,  _Quantum transport_ (Cambridge university press, 2009).
  * Likharev [1987] K. Likharev, Single-electron transistors: Electrostatic analogs of the dc SQUIDs, IEEE Trans. Mag. 23, 1142 (1987).
  * Fulton and Dolan [1987] T. A. Fulton and G. J. Dolan, Observation of single-electron charging effects in small tunnel junctions, Phys. Rev. Lett. 59, 109 (1987).
  * Nazarov [1999] Y. V. Nazarov, Coulomb blockade without tunnel junctions, Phys. Rev. Lett. 82, 1245 (1999).
  * Jezouin  _et al._ [2016] S. Jezouin, Z. Iftikhar, A. Anthore, F. D. Parmentier, U. Gennser, A. Cavanna, A. Ouerghi, I. P. Levkivskyi, E. Idrisov, E. V. Sukhorukov, L. I. Glazman, and F. Pierre, Controlling charge quantization with quantum fluctuations, Nature 536, 58 (2016).
  * Tinkham [1996] M. Tinkham,  _Introduction to superconductivity 2nd ed._ (McGraw-Hill international editions, 1996).
  * Blonder  _et al._ [1982] G. Blonder, M. Tinkham, and T. Klapwijk, Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion, Phys. Rev. B 25, 4515 (1982).
  * Pannetier and Courtois [2000] B. Pannetier and H. Courtois, Andreev reflection and proximity effect, J. Low Temp. Phys. 118, 599 (2000).
  * Fitzgerald  _et al._ [1998] R. J. Fitzgerald, S. L. Pohlen, and M. Tinkham, Observation of Andreev reflection in all-superconducting single-electron transistors, Phys. Rev. B 57, R11073 (1998).
  * Averin and Likharev [1986] D. V. Averin and K. K. Likharev, Coulomb blockade of single-electron tunneling, and coherent oscillations in small tunnel junctions, J. Low Temp. Phys. 62, 345 (1986).
  * Nakamura  _et al._ [1996] Y. Nakamura, C. D. Chen, and J. S. Tsai, Quantitative analysis of Josephson-quasiparticle current in superconducting single-electron transistors, Phys. Rev. B 53, 8234 (1996).
  * Pohlen  _et al._ [2000] S. L. Pohlen, R. J. Fitzgerald, and M. Tinkham, The Josephson–quasiparticle (jqp) current cycle in the superconducting single-electron transistor, Physica B 284, 1812 (2000).
  * Manninen  _et al._ [1997] A. J. Manninen, Y. A. Pashkin, A. N. Korotkov, and J. P. Pekola, Observation of thermally excited charge transport modes in a superconducting single-electron transistor, Europhys. Lett. 39, 305 (1997).
  * Lorenz  _et al._ [2018] T. Lorenz, S. Sprenger, and E. Scheer, Coulomb blockade and multiple Andreev reflection in a superconducting single-electron transistor, J. Low Temp. Phys. 191, 301 (2018).
  * Hekking  _et al._ [1995] F. W. J. Hekking, L. I. Glazman, and G. Schön, Andreev spectroscopy of Josephson coupling, Phys. Rev. B 51, 15312 (1995).
  * Fazio  _et al._ [2000] R. Fazio, G. Schön, and J. Siewert, Current–voltage characteristics of NSS transistors, Physica B 284-288, 1828 (2000).
  * Van Ruitenbeek  _et al._ [1996] J. van Ruitenbeek, A. Alvarez, I. Piñeyro, C. Grahmann, P. Joyez, M. Devoret, D. Esteve, and C. Urbina, Adjustable nanofabricated atomic size contacts, Rev. Sci. Instrum. 67, 108 (1996).
  * Scheer  _et al._ [1997] E. Scheer, P. Joyez, D. Esteve, C. Urbina, and M. H. Devoret, Conduction channel transmissions of atomic-size aluminum contacts, Phys. Rev. Lett. 78, 3535 (1997).
  * Agraıt  _et al._ [2003] N. Agraıt, A. L. Yeyati, and J. M. van Ruitenbeek, Quantum properties of atomic-sized conductors, Phys. Rep. 377, 81 (2003).
  * SM [2023] Supplemental material (2023), available online at….See supplemental material for further details.
  * Note [1] We show here data for positive $V_{\text{b}}$, but note that all maps were symmetric with respect to the bias polarity when choosing the adequate $n_{\rm g}$. For intermediate $\mathcal{T}_{\text{BJ}}\sim 0.5$ the peak slightly above $\Delta$ is hysteretic with respect to the bias sweep direction.
  * Nakamura  _et al._ [1995] Y. Nakamura, T. Sakamoto, and J. S. Tsai, Study of Josephson-quasiparticle cycles in superconducting single-electron transistors, JJAP 34, 4562 (1995).
  * Ambegaokar and Baratoff [1963] V. Ambegaokar and A. Baratoff, Tunneling between superconductors, Phys. Rev. Lett. 10, 486 (1963).
  * Joyez  _et al._ [1994] P. Joyez, P. Lafarge, A. Filipe, D. Esteve, and M. Devoret, Observation of parity-induced suppression of Josephson tunneling in the superconducting single electron transistor, Phys. Rev. Lett. 72, 2458 (1994).
  * Cleuziou  _et al._ [2006] J.-P. Cleuziou, W. Wernsdorfer, V. Bouchiat, T. Ondarçuhu, and M. Monthioux, Carbon nanotube superconducting quantum interference device, Nature nanotechnology 1, 53 (2006).
  * van den Brink  _et al._ [1991] A. M. van den Brink, G. Schön, and L. Geerligs, Combined single-electron and coherent-Cooper-pair tunneling in voltage-biased Josephson junctions, Phys. Rev. Lett. 67, 3030 (1991).
  * Averin and Aleshkin [1989] D. V. Averin and V. Y. Aleshkin, Resonance tunneling of Cooper pairs in a system of two small Josephson junctions, JETP Lett, 50, 367 (1989).
  * Cuevas and Belzig [2003] J. C. Cuevas and W. Belzig, Full counting statistics of multiple Andreev reflections, Phys. Rev. Lett. 91, 187001 (2003).
  * Note [2] For describing the cotunneling processes at intermediate coupling we will also take a small AR contribution through the TB junction into account.
  * Pekola  _et al._ [2010] J. P. Pekola, V. F. Maisi, S. Kafanov, N. Chekurov, A. Kemppinen, Y. A. Pashkin, O.-P. Saira, M. Möttönen, and J. S. Tsai, Environment-assisted tunneling as an origin of the Dynes density of states, Phys. Rev. Lett. 105, 026803 (2010).
  * Maisi  _et al._ [2011] V. F. Maisi, O. P. Saira, Y. A. Pashkin, J. S. Tsai, D. V. Averin, and J. P. Pekola, Real-time observation of discrete Andreev tunneling events, Phys. Rev. Lett. 106, 217003 (2011).
  * Note [3] We note that Al single-atom contacts with $\mathcal{T}\gtrsim 0.8$ are expected to accommodate more than one channel [19]. For this contact with $\mathcal{T}=0.4$, it was sufficient to consider one single channel with $\tau=0.4$. Considering more than one channel would result in a smaller current, and hence, would not improve the agreement between experiment and simulation giving a worse agreement.
  * Joyez [1995] P. Joyez,  _Le transistor a une paire de Cooper: Un système quantique macroscopique_ , Ph.D. thesis, Université Paris 6, CEA-Saclay (1995).
  * Siewert and Schön [1996] J. Siewert and G. Schön, Charge transport in voltage-biased superconducting SET transistors, Phys. Rev. B 54, 7421 (1996).
  * Siddiqi [2021] I. Siddiqi, Engineering high-coherence superconducting qubits, Nature Reviews Materials 6, 875–891 (2021).