# Quantum catalysis in cavity QED

A. de Oliveira Junior  Faculty of Physics, Astronomy and Applied Computer
Science, Jagiellonian University, 30-348 Kraków, Poland.  Martí Perarnau-
Llobet  Department of Applied Physics, University of Geneva, 1211 Geneva,
Switzerland  Nicolas Brunner  Department of Applied Physics, University of
Geneva, 1211 Geneva, Switzerland  Patryk Lipka-Bartosik  Department of Applied
Physics, University of Geneva, 1211 Geneva, Switzerland

###### Abstract

Catalysis plays a key role in many scientific areas, most notably in chemistry
and biology. Here we present a catalytic process in a paradigmatic quantum
optics setup, namely the Jaynes-Cummings model, where an atom interacts with
an optical cavity. The atom plays the role of the catalyst, and allows for the
deterministic generation of non-classical light in the cavity. Considering a
cavity prepared in a “classical” coherent state, and choosing appropriately
the atomic state and the interaction time, we obtain an evolution with the
following properties. First, the state of the cavity has been modified, and
now features non-classicality, as witnessed by sub-Poissonian statistics or
Wigner negativity. Second, the process is catalytic, in the sense that the
atom is deterministically returned to its initial state exactly, and could
then in principle be re-used multiple times. We investigate the mechanism of
this catalytic process, in particular highlighting the key role of
correlations and quantum coherence.

##  I Introduction

The effect of catalysis involves using an auxiliary system (a catalyst) to
enable a process that would either not occur spontaneously or would occur very
slowly. Catalysis manifests across a variety of fields (see e.g. Chorkendorff
and Niemantsverdriet (2017)), including biological processes activated by
enzymes, the speed-up of chemical reactions, and the synthesis of
nanomaterials.

More recently, the phenomenon of catalysis has also become relevant in the
context of quantum information; see recent reviews Datta  _et al._ (2022a);
Lipka-Bartosik  _et al._ (2023a). First examples focused on entanglement
manipulation Jonathan and Plenio (1999); van Dam and Hayden (2003); Turgut
(2007); Daftuar and Klimesh (2001); Sun  _et al._ (2005); Feng  _et al._
(2005); Lipka-Bartosik and Skrzypczyk (2021a); Kondra  _et al._ (2021); Datta
_et al._ (2022b), and then spread to other areas, including quantum
thermodynamics Brandão  _et al._ (2015); Ng  _et al._ (2015); Wilming and
Gallego (2017); Müller (2018); Lipka-Bartosik and Skrzypczyk (2021b);
Shiraishi and Sagawa (2021); Gallego  _et al._ (2016); Boes  _et al._ (2020);
Henao and Uzdin (2021, 2023); Son and Ng (2022, 2023); Lipka-Bartosik  _et
al._ (2023b); Czartowski  _et al._ (2023), coherence theory Åberg (2014);
Vaccaro  _et al._ (2018); Lostaglio and Müller (2019); Takagi and Shiraishi
(2022); Char  _et al._ (2023); van Luijk  _et al._ (2023) and others Marvian
and Spekkens (2019); Wilming (2021, 2022); Rubboli and Tomamichel (2022); Lie
and Jeong (2021); Boes  _et al._ (2019). These results are typically
formulated within the framework of quantum resource theories Chitambar and
Gour (2019). This abstract approach is particularly useful for characterizing
the fundamental limits of manipulating quantum resources, including scenarios
involving catalytic systems of arbitrary complexity.

An interesting direction is whether quantum catalysis is also relevant and
useful in a more practical context, potentially even in experiments. Here we
investigate quantum catalysis in a paradigmatic setup of quantum optics,
namely the Jaynes-Cummings model Jaynes and Cummings (1963); Larson and
Mavrogordatos (2021); Greentree  _et al._ (2013), where a two-level atom
interacts with a single-mode optical cavity. We uncover a catalytic process
enabling the generation of a non-classical state of light in the cavity, using
the atom as a catalyst. Specifically, we consider the cavity to be initially
prepared in a “classical” coherent state, and uncorrelated to the atom. By
carefully setting the initial state of the atom and the interaction time, we
obtain a final state such that (i) the atom is back in its initial state
exactly, and (ii) the state of the cavity is now non-classical, i.e. featuring
Wigner negativity or sub-Poisonian statistics. Hence, non-classicalilty of the
cavity has been generated without perturbing the state of the atom (see Fig.
1). The process is catalytic and the atom could be re-used, for example by
coupling it to another cavity.

Figure 1: Quantum catalysis in the Jaynes-Cummings model. (a) An atom (the
catalyst C) interacts with a single-mode optical cavity (the system S),
initially prepared in a “classical” coherent state. (b) We consider the
evolution $U$ over a well-chosen time interval (from $t=0$ to $t=\tau$) such
that (i) the final state of the cavity is non-classical, and (ii) the atom is
returned to its intial state exactly. Hence, non-classicality has been
generated via catalysis.

We investigate the mechanism of this catalytic process, and identify two
crucial ingredients. First, the final state of the atom and cavity must
feature correlations. Second, the evolution of the state of the catalyst must
involve quantum coherence (i.e. superpositions of the energy basis states).
The latter is an interesting aspect, as typical instances of quantum catalysis
in resource theories involve only diagonal states (i.e. without coherence), so
that they can be understood as a stochastic process involving the probability
distributions of the system and catalysis. Instead here, the system and
catalyst experience a genuinely quantum evolution. This is a novel instance of
the effect of _coherent quantum catalysis_ , recently investigated in quantum
thermodynamics Lipka-Bartosik  _et al._ (2023b).

Before proceeding, it is worth discussing previous works in quantum optics
that relate to the concept of catalysis. Notably, the pioneering proposal for
quantum computing in ion traps Cirac and Zoller (1995) (see also Phoenix and
Barnett (1993); Hagley  _et al._ (1997)) considers two spin qubits that become
entangled via an interaction with a cavity that can be considered catalytic
Lipka-Bartosik  _et al._ (2023a). Another relevant direction is that of
“multi-photon catalysis” (see e.g. Lvovsky and Mlynek (2002); Bartley  _et
al._ (2012); Hu  _et al._ (2016, 2017)) which is a heralded catalytic process,
where the catalyst is returned only with some probability. In contrast, our
catalytic protocol is deterministic.

##  II Catalysis in the Jaynes-Cummings model

Let us first define the concept of quantum catalysis. Consider a quantum setup
comprising a system $(\textsf{S})$ initially prepared in a state
$\rho_{\textsf{S}}$ and a catalyst $(\textsf{C})$ in an initial state
$\chi_{\textsf{C}}$. The total system SC is assumed to be closed and evolves
via an energy-conserving process for some time $\tau$. This evolution is
represented by the unitary $U=\exp({-iH_{\textsf{SC}}\tau})$, where
$H_{\textsf{SC}}$ denotes the Hamiltonian. Consequently, the final state of
the total system is given by
$\sigma_{\textsf{SC}}:=U(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})U^{\dagger}$.

The evolution is said to be catalytic when the catalyst is returned in exactly
the same state as it was initially prepared. Formally, we demand that

$\sigma_{\textsf{C}}:=\Tr_{\textsf{S}}[U(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})U^{\dagger}]=\chi_{\textsf{C}},$
(1)

which we refer to as the _catalytic constraint_. Satisfying this constraint
typically requires to carefully choose the initial states of the system
$\rho_{\textsf{S}}$ and the catalyst $\chi_{\textsf{C}}$, as well as the
interaction time $\tau$.

The main goal of a catalytic evolution is to induce an interesting local
dynamics on the system S, i.e.

$\displaystyle\rho_{\textsf{S}}\rightarrow\sigma_{\textsf{S}}:=\Tr_{\textsf{C}}[U(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})U^{\dagger}],$
(2)

while leaving the state of the catalyst unchanged. Notably, it is possible to
induce an evolution on S [as in Eq. (2)] that would not be possible without
the presence of the catalyst.

In this work, we discuss the phenomenon of catalysis in the Jaynes-Cummings
(JC) model, describing the interaction between a single-mode optical cavity
and a two-level atom Jaynes and Cummings (1963) (see Fig. 1a). Here we choose
the cavity to represent the system S, while the atom will play the role of the
catalyst C. The cavity is characterized by the bosonic creation and
annihilation operators $a^{\dagger}$ and $a$ with the photon number operator
$n_{\textsf{S}}:=a^{\dagger}a$. The atom has energy levels $\ket{g}$ and
$\ket{e}$ and its energy is captured by the operator
$\sigma_{z}=\outerproduct{e}{e}-\outerproduct{g}{g}$. We work in the resonant
regime, where the atom and cavity have the same frequency $\omega$. The
evolution is governed by the JC Hamiltonian, which in the rotating-wave
approximation reads

$\displaystyle H_{\textsf{SC}}=\omega
a^{\dagger}a+\frac{\omega}{2}\sigma_{z}+g\left(\sigma_{+}a+\sigma_{-}a^{\dagger}\right),$
(3)

where $g$ is the coupling constant and $\sigma_{+}=\outerproduct{e}{g}$ and
$\sigma_{-}=\outerproduct{g}{e}$. Note that, as we focus on the resonant
regime, the evolution specified by Eq. (3) is energy preserving.

The task we would like to achieve in a catalytic way is the generation of non-
classical light in the cavity. We consider an initial state of the cavity that
is classical, namely a coherent state

$\ket{\alpha}=e^{-|\alpha|^{2}/2}\sum_{n=0}^{\infty}\frac{\alpha^{n}}{n!}\ket{n}.$
(4)

This state has Poissonian statistics with the mean number of photons $\langle
n_{\textsf{S}}\rangle_{\outerproduct{\alpha}{\alpha}}=|\alpha|^{2}$.

Figure 2: First illustrative example. Catalytic process for generating non-
classicality in the cavity, as captured by the second-order auto-correlation
function. Its time evolution, $g^{(2)}(t):=g^{(2)}(\sigma_{\textsf{S}}(t))$,
is shown in the top panel, while the bottom panel shows the modification of
the atomic state, via the distance
$\Delta(t):=\norm{{\chi_{\textsf{C}}-\sigma_{\textsf{C}}(t)}}_{1}$ with
respect to the initial state. The orange stars indicate the final time
($\tau\approx 40$) for which the evolution is catalytic. The atom returns to
its intial state ($\Delta(\tau)=0$), while non-classicality has been activated
$g^{(2)}(\tau)\approx 0.5<g^{(2)}(0)=1$. The value of the parameters are:
$\alpha=1/\sqrt{2}$, $\omega=2\pi$, $g=\pi$.

Our goal is to find an initial state of the atom $\chi_{\textsf{C}}$ and an
interaction time $\tau$ such that the evolution is catalytic [i.e. satisfying
Eq. (1)], while at the same time leading to a final state of the cavity
$\sigma_{\textsf{S}}$ that is non-classical. Below we present two illustrative
examples of such a process, using two complementary figures of merit to
witness non-classicality. Note that the detailed analysis of these examples
will be discussed in the next section.

Figure 3: Second illustrative example. Catalytic generation of Wigner
negativity. Panel $(a)$ shows the time evolution of the Wigner logarithmic
negativity $\mathsf{W}$, as well as the Wigner function at initial and final
time $\tau$. Non-classicality is clearly generated, while the process is
catalytic. Panel $(b)$ shows the trajectory of the atomic state, i.e. the
catalyst C, on the Bloch sphere. The initial and final state coincide (red
dot). The values of the parameters are: $\alpha=1/\sqrt{2}$, $\omega=2\pi$ and
$g=\pi$.

For our first example, we quantify non-classicality via the second-order auto
correlation function of the final state of the cavity Glauber (1963) , i.e.

$\displaystyle g^{(2)}(\sigma_{\textsf{S}})=\frac{\langle
n_{\textsf{S}}^{2}\rangle_{\sigma}-\langle
n_{\textsf{S}}\rangle_{\sigma}}{\langle n_{\textsf{S}}\rangle_{\sigma}^{2}},$
(5)

where $\langle
n_{\textsf{S}}\rangle_{\rho}:=\Tr[\rho_{\textsf{S}}n_{\textsf{S}}]$. Our goal
is to obtain $g^{(2)}(\sigma_{\textsf{S}})<1$, which certifies a non-classical
state of light.

In Fig. 2, we show such a catalytic activation of non-classicality. The upper
plot depicts the time evolution of $g^{(2)}$. Catalysis occurs at time
$\tau\approx 40$, which is represented by the orange stars. Starting from the
initial state of the cavity
$\rho_{\textsf{S}}=\outerproduct{\alpha}{\alpha}_{\textsf{S}}$ (with
$\alpha=1/\sqrt{2}$), which has $g^{(2)}(\rho_{\textsf{S}})=1$ as any coherent
state, we obtain a final state $\sigma_{\textsf{S}}$ for which
$g^{(2)}(\sigma_{\textsf{S}})\approx 0.5$, hence non-classical. In parallel,
we monitor the time evolution of the atomic state by calculating the trace
distance to its initial state, i.e.

$\displaystyle\Delta(t):=\norm{{\chi_{\textsf{C}}-\sigma_{\textsf{C}}(t)}}_{1}.$
(6)

We see that at the final time $\tau\approx 40$ we get $\Delta(\tau)=0$. Hence
the atom returned to its initial state exactly, as required for catalysis.

As a second example, we demonstrate how Wigner negativity can be catalytically
generated. We focus on the Wigner logarithmic negativity Veitch  _et al._
(2014) (WLN), defined as

$\displaystyle\mathsf{W}\left(\rho\right):=\log\left(\int\\!\mathrm{d}x\,\mathrm{d}p\,\left|W_{\rho}\left(x,p\right)\right|\right),$
(7)

where $W_{\rho}\left(x,p\right)$ is the Wigner function

$\displaystyle W_{\rho}(x,p)=\frac{1}{\pi}\int e^{2ipx^{\prime}}\langle
x-x^{\prime}|\rho|x+x^{\prime}\rangle\,\text{d}x^{\prime}.$ (8)

Starting from a coherent state
$\rho_{\textsf{S}}=\outerproduct{\alpha}{\alpha}_{\textsf{S}}$ which has a
positive Wigner function, hence $\mathsf{W}(\rho_{\textsf{S}})=0$, we aim to
obtain a final state of the cavity $\sigma_{\textsf{S}}$ with
$\mathsf{W}(\sigma_{\textsf{S}})>0$, therefore certifying its non-
classicality. In Fig. 3, we present an example of such an evolution. First, we
plot WLN as a function of time $t$. At the final time $\tau\approx 5$, we
obtain a state $\sigma_{\textsf{S}}$ that has
$\mathsf{W}(\sigma_{\textsf{S}})\approx 0.1$, and we plot its Wigner function.
To verify the catalytic nature of the evolution, we display the evolution of
the atomic state via its trajectory in the Bloch sphere. Crucially, the
trajectory is closed, as the initial and final state of the atom exactly
coincide (red dot).

Before moving to the next section, we notice that the two examples we analysed
are complementary. The first one shows the generation of non-classicality as
witnessed by the $g^{(2)}$ function, while the final state
$\sigma_{\textsf{S}}$ has a positive Wigner function. In the second example,
the opposite happens. Namely, we observe that the final state
$\sigma_{\textsf{S}}$ has a negative Wigner function, even though
$g^{(2)}(\sigma)>1$.

##  III Mechanism of catalysis

Here we start by providing a more intuitive understanding of quantum
catalysis, by identifying a mechanism which allows for the activation of non-
classicality. In turn, this allows us to characterize analytically the
catalytic regime.

In the Jaynes-Cummings model, the energies of the cavity S and the atom C are
specified by local number operators $n_{\textsf{S}}$ and $n_{\textsf{C}}$
respectively. Hence, the total energy of both systems is proportional to the
number of excitations, and described by a joint operator
$n_{\textsf{SC}}:=n_{\textsf{S}}+n_{\textsf{C}}$. Since the JC evolution
$U(t)$ that takes $\rho_{\textsf{S}}\otimes\chi_{\textsf{C}}$ into
$\sigma_{\textsf{SC}}$ conserves the total energy, we have that
$[U(t),n_{\textsf{S}}+n_{\textsf{C}}]=0$ for all $t$. When the evolution is
catalytic, all moments of $n_{\textsf{C}}$ must remain unchanged, in
particular $\langle n_{\textsf{C}}\rangle_{\chi}$ = $\langle
n_{\textsf{C}}\rangle_{\sigma}$. Consequently, the first moment of
$n_{\textsf{S}}$ is also preserved, that is $\langle
n_{\textsf{S}}\rangle_{\rho}=\langle n_{\textsf{S}}\rangle_{\sigma}$.
Importantly, this is not the case for higher moments of $n_{\textsf{S}}$.
Particularly, in Appendix B we show that the second moment satisfies

$\displaystyle\langle n^{2}_{\textsf{S}}\rangle_{\sigma}=\langle
n^{2}_{\textsf{S}}\rangle_{\rho}+2\Bigg{(}\langle
n_{\textsf{S}}\rangle_{\sigma}\langle n_{\textsf{C}}\rangle_{\sigma}-\langle
n_{\textsf{S}}\otimes n_{\textsf{C}}\rangle_{\sigma}\Bigg{)},$ (9)

Hence, the second moment in the final state of the cavity, $\langle
n^{2}_{\textsf{S}}\rangle_{\sigma}$, can become smaller (or larger) than the
second moment in the initial state $\langle n^{2}_{\textsf{S}}\rangle_{\rho}$.
This means that using a catalyst allows for modifying the distribution of the
local observable $n_{\textsf{S}}$ of the system S: while its average must
remain the same, the higher moments can change. Importantly, this can only
happen if the system becomes correlated with the catalyst, i.e. $\langle
n^{2}_{\textsf{S}}\rangle_{\sigma}\neq\langle
n^{2}_{\textsf{S}}\rangle_{\rho}$ only if
$\sigma\neq\sigma_{\textsf{S}}\otimes\sigma_{\textsf{C}}$ [as seen from Eq.
(9)]. Thus, these correlations are essential for observing quantum catalysis.
Finally, we note that the above analysis also applies beyond the JC model, see
Appendix B-1a for a generalization of Eq. (9) to arbitrary observables and
moments. This can help to discover and understand new instances of quantum
catalysis.

The above analysis will also serve as a basis for the characterisation of the
parameter regime leading to catalysis. In particular we derive a necessary and
sufficient condition for satisfying the catalytic constraint of Eq. (1).
Moreover, we obtain an analytic expression for the auto-correlation function
$g^{2}$, which is based on the second moment $\langle
n^{2}_{\textsf{S}}\rangle_{\sigma}$.

We consider an arbitrary initial state of the atom

$\displaystyle\chi_{\textsf{C}}=q\left|g\right\rangle\\!\\!\left\langle
g\right|+r\left|g\right\rangle\\!\\!\left\langle
e\right|+r^{*}\left|e\right\rangle\\!\\!\left\langle
g\right|+[1-q]\left|e\right\rangle\\!\\!\left\langle e\right|,$ (10)

as well as a general initial state of the cavity
$\rho_{\textsf{S}}=\sum_{n,m}^{\infty}p_{n,m}\outerproduct{n}{m}$ with
$p_{n}:=p_{n,n}$. Combining Eq. (9) with the fact that mean energy of the
cavity is conserved, i.e. $\langle n_{\textsf{S}}\rangle_{\sigma}=\langle
n_{\textsf{S}}\rangle_{\rho}$, we get

$\displaystyle\\!g^{(2)}(\sigma_{\textsf{S}})=g^{(2)}(\rho_{\textsf{S}})-\frac{2}{\langle
n_{\textsf{S}}\rangle_{\rho}^{2}}\\!\left[\langle n_{\textsf{S}}\otimes
n_{\textsf{C}}\rangle_{\sigma}\\!-(1-q)\langle
n_{\textsf{S}}\rangle_{\rho}\right]\\!,\\!$ (11)

where we have

$\displaystyle\\!\\!\\!\langle n_{\textsf{S}}\otimes
n_{\textsf{C}}\rangle_{\sigma}\\!=\\!\sum_{n=0}^{\infty}n\left[(1-q)p_{n}c_{n}^{2}+y_{n}+qp_{n+1}s_{n}^{2}\right],$
(12)

with $s_{n}:=\sin(gt\sqrt{n+1})$, $c_{n}:=\cos(gt\sqrt{n+1})$, and
$y_{n}:=2\operatorname{Im}\quantity[rp_{n+1,n}]s_{n}c_{n}$. For details see
Appendix A-3a.

In order to satisfy the catalytic constraint, we obtain a set of equations for
the components of the atomic state. Decomposing the diagonal term as
$q=q_{\text{inc}}+q_{\text{coh}}$, we get

$\displaystyle
q_{\text{inc}}=\frac{1}{Q}\sum_{n=0}^{\infty}p_{n}s_{n}^{2},\qquad
q_{\text{coh}}=\frac{1}{Q}\sum_{n=0}^{\infty}y_{n},$ (13)

with $Q:=\sum_{n=0}^{\infty}(p_{n}+p_{n+1})s_{n}^{2}$. Interestingly,
$q_{\text{inc}}$ is specified by the occupations of $\rho_{\textsf{S}}$, while
$q_{\text{coh}}$ depends on its coherence in the Fock basis. Moreover, the
off-diagonal term $r$ satisfies

$r=\frac{i(a_{3}a^{*}_{4}+a^{*}_{1}a_{4})}{|a_{1}|^{2}-|a_{3}|^{2}}-\frac{i(a_{3}a^{*}_{2}+a^{*}_{1}a_{2})}{|a_{1}|^{2}-|a_{3}|^{2}}q,$
(14)

with $a_{i}$ being auxiliary functions defined as

$\displaystyle\begin{split}a_{1}&=\sum_{n=0}^{\infty}p_{n,n}c_{n-1}c_{n}-e^{-i\omega\tau},\qquad
a_{3}=\sum_{n=0}^{\infty}p_{n,n+2}s_{n}s_{n+1},\\\
a_{2}&=\sum_{n=0}^{\infty}p_{n,n+1}s_{n}\quantity[c_{n-1}+c_{n+1}],\hskip
10.0pta_{4}=\sum_{n=0}^{\infty}p_{n,n+1}s_{n}c_{n+1}.\end{split}$ (15)

For a detailed derivation of Eqs. (15) see Appendix A-1.

Importantly, Eqs (13) and (14) are necessary and sufficient for ensuring the
catalytic constraint of Eq. (1). In combination with Eq. (11), we can now
characterize analytically the effect of non-classicality activation for the
$g^{2}$ function in the catalytic regime. This is done for finding the
parameters (initial states and interaction time) for the first illustrative
example of a catalytic process (see Fig. 2). For the second example, we use
again such analytic results to ensure the validity of the catalytic
constraint, while the Wigner functions are computed numerically.

##  IV How general is catalysis?

Figure 4: Which states lead to catalysis? Panel (a) shows the minimal value of
$g^{2}$ obtained as a function of the amplitude $|\alpha|$ of the initial
coherent state of the cavity
$\rho_{\textsf{S}}=\outerproduct{\alpha}{\alpha}_{\textsf{S}}$. The orange
star corresponds to our first illustrative example. Panel (b) displays the
atomic states (in the $y-z$ plane of the Bloch sphere) that satisfy the
catalytic constraint and generate non-classical states, for an intial coherent
state $\alpha=1/\sqrt{2}$. The colour represents different values of
$g^{(2)}<1$. We impose a limit on the interaction time $\tau\leq 100$ and take
$10^{6}$ samples. In both panels, parameters are $\omega=2\pi$, $g=\pi$.

An interesting problem is to understand how typical the effect of catalysis
is. Here we discuss different aspects of this question.

To begin with, note that it is not obvious a priori whether the catalytic
constraint of Eq. (1) can be satisfied. However, due to the quantum version of
the Perron-Frobenius theorem, every quantum channel has at least one positive
semi-definite fixed point Fannes  _et al._ (1992). Now, for a fixed input
state $\rho_{\textsf{S}}$ and a fixed interaction time $\tau$, the state of
the atom C evolves according to an effective quantum channel
$\chi_{\textsf{C}}\rightarrow\Tr_{\textsf{S}}[U(\tau)(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})U^{\dagger}(\tau)]$.
Consequently, there always exists an initial state $\chi_{\textsf{C}}$ which
is left unchanged by this channel, hence providing at least one solution to
Eq. (1).

Next, one might wonder how often a catalytic evolution leads to a non-
classical state of the cavity. In particular, when preparing the cavity in a
coherent state $\ket{\alpha}$, does there always exist a state of the catalyst
$\chi_{\textsf{C}}$ which allows to generate non-classicality? To address this
question, we investigate the minimum value of $g^{(2)}$ of the final state
$\sigma_{\textsf{S}}$, as a function of $|\alpha|$ [see Fig. 4(a)]. This is
done by combining Eqs. (11-12) and (13-14), and imposing a bound on the final
time, i.e $g\tau\leq 100$. For $\alpha\in(0,2]$, we observe that
$g^{(2)}(\sigma_{\textsf{S}})<1$, indicating that non-classicality generation
via catalysis is generic here. Note that when the initial coherent state has
low energy, the final state of the cavity is close to the intial one, but with
a slightly reduced variance, leading to a value of $g^{(2)}$ approaching zero.

Let us now ask the converse question, i.e. whether every atomic state can lead
to a catalytic evolution. A first observation is that pure states cannot act
as useful catalysts in general, and in particular cannot catalytically
generate non-classicality. Indeed, a key ingredient for catalysis is the fact
that the system and catalyst become correlated [see Eq. (9)], which is
impossible when the state of the catalyst is pure.

To further explore this question, we also investigated which states of the
atom can catalytically generate sub-Poissonian statistics (i.e. $g^{(2)}<1$),
given an initial coherent state of the cavity and a limited interaction time
$g\tau\leq 100$. In Fig. 4(b), we display an example of such a set of
catalytic states. Interestingly, this set appears to contain states that are
almost pure. The structure of the catalytic set is further discussed in the
Appendix A-1, in particular we observe a strong dependence with the intial
state of the cavity.

Finally, we also investigated the case where the initial state of the cavity
is not a coherent one. Interestingly, we find that catalysis can boost non-
classicality. Specifically, starting with an initial state $\rho_{\textsf{S}}$
that has some level of non-classicality, one can get a final state
$\sigma_{\textsf{S}}$ that features more non-classicality, as quantified by
the auto-correlation function, i.e.
$g^{(2)}(\sigma_{\textsf{S}})<g^{(2)}(\rho_{\textsf{S}})$ (see Appendix C-2
for details). Conversely, there also exist initial states of the cavity (e.g.
Fock states) for which non-classicality cannot be increased catalytically (see
Appendix C-1). This is intuitive, as Fock states are the most non-classical
(i.e. with minimal value of $g^{(2)}$) among all states with a fixed average
energy. As a catalytic protocol must preserve the mean energy, the value of
$g^{(2)}$ cannot be increased.

##  V Discussion

We presented a catalytic process for generating non-classical states of light
in an optical cavity. Our work shows that quantum catalysis, a concept so far
explored in the abstract framework of resource theories, is relevant in a
practical context. Furthermore, our protocol could potentially be implemented
in state-of-the-art experimental setups Frisk Kockum  _et al._ (2019); Periwal
_et al._ (2021); de Buy Wenniger  _et al._ (2022), e.g. in cavity QED
Ballester  _et al._ (2012); Blais  _et al._ (2021) or trapped ions Leibfried
_et al._ (2003); Lv  _et al._ (2018). Beyond proof-of-principle experiments,
it would also be interesting to investigate whether such a catalytic protocol
offers a practical advantage. Indeed, the key point of catalysis is that the
catalyst (here the atom) is returned exactly in the same state as it was
initially prepared. Hence the same atom could in principle be used repeatedly
for activating non-classicality in different cavities, or in the same cavity,
but at different times. In the future, it would be interesting to uncover
further instances of quantum catalysis in realistic setups, e.g. exploring
various platforms and models in quantum optics, as well as for applications in
quantum information and metrology.

Another relevant aspect of our work is the fact that the catalyst (atom) is,
in general, a coherent quantum system. Specifically, the activation of non-
classicality requires a state of the atom that, during its evolution, features
coherence with respect to the energy basis. In this sense, it is an instance
of coherent quantum catalysis Lipka-Bartosik  _et al._ (2023b). This contrasts
with most previous examples of quantum catalysis, and a deeper understanding
of the role of coherence in catalysis is an exciting future direction.

##  VI Acknowledgements

We thank Henrik Wilming and Nelly H.Y. Ng for valuable discussions. AOJ
acknowledge financial support from the Foundation for Polish Science through
the TEAM-NET project (contract no. POIR.04.04.00-00-17C1/18-00); MPL, NB and
PLB acknowledge the Swiss National Science Foundation for financial support
through the Ambizione grant PZ00P2-186067 and the NCCR SwissMAP.

## References

  * Chorkendorff and Niemantsverdriet (2017) I. Chorkendorff and J. W. Niemantsverdriet,  _Concepts of modern catalysis and kinetics_ (John Wiley & Sons, 2017).
  * Datta  _et al._ (2022a) C. Datta, T. V. Kondra, M. Miller, and A. Streltsov, arXiv:2209.15213 (2022a).
  * Lipka-Bartosik  _et al._ (2023a) P. Lipka-Bartosik, H. Wilming, and N. H. Y. Ng, Catalysis in quantum information theory, _in preparation_ (2023a).
  * Jonathan and Plenio (1999) D. Jonathan and M. B. Plenio, Phys. Rev. Lett. 83, 3566–3569 (1999).
  * van Dam and Hayden (2003) W. van Dam and P. Hayden, Phys. Rev. A 67, 060302 (2003).
  * Turgut (2007) S. Turgut, J. Phys. A Math. Theor. 40, 12185 (2007).
  * Daftuar and Klimesh (2001) S. Daftuar and M. Klimesh, Phys. Rev. A. 64, 042314 (2001).
  * Sun  _et al._ (2005) X. Sun, R. Duan, and M. Ying, IEEE Trans. Inf. Theory 51, 75 (2005).
  * Feng  _et al._ (2005) Y. Feng, R. Duan, and M. Ying, IEEE Trans. Inf. Theory 51, 1090 (2005).
  * Lipka-Bartosik and Skrzypczyk (2021a) P. Lipka-Bartosik and P. Skrzypczyk, Phys. Rev. Lett. 127, 080502 (2021a).
  * Kondra  _et al._ (2021) T. V. Kondra, C. Datta, and A. Streltsov, Phys. Rev. Lett. 127, 150503 (2021).
  * Datta  _et al._ (2022b) C. Datta, T. V. Kondra, M. Miller, and A. Streltsov, arXiv:2207.05694 (2022b).
  * Brandão  _et al._ (2015) F. Brandão, M. Horodecki, N. Ng, J. Oppenheim, and S. Wehner, PNAS 112, 3275–3279 (2015).
  * Ng  _et al._ (2015) N. H. Y. Ng, L. Mančinska, C. Cirstoiu, J. Eisert, and S. Wehner, New J. Phys. 17, 085004 (2015).
  * Wilming and Gallego (2017) H. Wilming and R. Gallego, Phys. Rev. X 7, 041033 (2017).
  * Müller (2018) M. P. Müller, Phys. Rev. X 8, 041051 (2018).
  * Lipka-Bartosik and Skrzypczyk (2021b) P. Lipka-Bartosik and P. Skrzypczyk, Phys. Rev. X 11, 011061 (2021b).
  * Shiraishi and Sagawa (2021) N. Shiraishi and T. Sagawa, Phys. Rev. Lett. 126, 150502 (2021).
  * Gallego  _et al._ (2016) R. Gallego, J. Eisert, and H. Wilming, New J. Phys. 18, 103017 (2016).
  * Boes  _et al._ (2020) P. Boes, R. Gallego, N. H. Ng, J. Eisert, and H. Wilming, Quantum 4, 231 (2020).
  * Henao and Uzdin (2021) I. Henao and R. Uzdin, Quantum 5, 547 (2021).
  * Henao and Uzdin (2023) I. Henao and R. Uzdin, Phys. Rev. Lett. 130, 020403 (2023).
  * Son and Ng (2022) J. Son and N. H. Y. Ng, arXiv:2209.15213 (2022).
  * Son and Ng (2023) J. Son and N. H. Y. Ng, arXiv:2303.13020 (2023).
  * Lipka-Bartosik  _et al._ (2023b) P. Lipka-Bartosik, M. Perarnau-Llobet, and N. Brunner, Phys. Rev. Lett. 130, 040401 (2023b).
  * Czartowski  _et al._ (2023) J. Czartowski, A. de Oliveira Junior, and K. Korzekwa, arXiv:2303.12840 (2023).
  * Åberg (2014) J. Åberg, Phys. Rev. Lett. 113, 150402 (2014).
  * Vaccaro  _et al._ (2018) J. A. Vaccaro, S. Croke, and S. M. Barnett, J. Phys. A Math. Theor. 51, 414008 (2018).
  * Lostaglio and Müller (2019) M. Lostaglio and M. P. Müller, Phys. Rev. Lett. 123, 020403 (2019).
  * Takagi and Shiraishi (2022) R. Takagi and N. Shiraishi, Phys. Rev. Lett. 128, 240501 (2022).
  * Char  _et al._ (2023) P. Char, D. Chakraborty, A. Bhar, I. Chattopadhyay, and D. Sarkar, Phys. Rev. A. 107, 012404 (2023).
  * van Luijk  _et al._ (2023) L. van Luijk, R. F. Werner, and H. Wilming, arXiv:2301.09877 (2023).
  * Marvian and Spekkens (2019) I. Marvian and R. W. Spekkens, Phys. Rev. Lett. 123, 020404 (2019).
  * Wilming (2021) H. Wilming, Phys. Rev. Lett. 127, 260402 (2021).
  * Wilming (2022) H. Wilming, Quantum 6, 858 (2022).
  * Rubboli and Tomamichel (2022) R. Rubboli and M. Tomamichel, Phys. Rev. Lett. 129, 120506 (2022).
  * Lie and Jeong (2021) S. H. Lie and H. Jeong, Phys. Rev. Res. 3, 043089 (2021).
  * Boes  _et al._ (2019) P. Boes, J. Eisert, R. Gallego, M. P. Müller, and H. Wilming, Phys. Rev. Lett. 122, 210402 (2019).
  * Chitambar and Gour (2019) E. Chitambar and G. Gour, Rev. Mod. Phys. 91, 025001 (2019).
  * Jaynes and Cummings (1963) E. T. Jaynes and F. W. Cummings, Proc. IEEE 51, 89 (1963).
  * Larson and Mavrogordatos (2021) J. Larson and T. Mavrogordatos,  _The Jaynes-Cummings Model and Its Descendants_ (IOP Publishing, 2021).
  * Greentree  _et al._ (2013) A. D. Greentree, J. Koch, and J. Larson, J. Phys. B 46, 220201 (2013).
  * Cirac and Zoller (1995) J. I. Cirac and P. Zoller, Phys. Rev. Lett. 74, 4091 (1995).
  * Phoenix and Barnett (1993) S. J. Phoenix and S. M. Barnett, Journal of Modern Optics 40, 979 (1993), https://doi.org/10.1080/09500349314551051 .
  * Hagley  _et al._ (1997) E. Hagley, X. Maître, G. Nogues, C. Wunderlich, M. Brune, J. M. Raimond, and S. Haroche, Phys. Rev. Lett. 79, 1 (1997).
  * Lvovsky and Mlynek (2002) A. I. Lvovsky and J. Mlynek, Phys. Rev. Lett. 88, 250401 (2002).
  * Bartley  _et al._ (2012) T. J. Bartley, G. Donati, J. B. Spring, X.-M. Jin, M. Barbieri, A. Datta, B. J. Smith, and I. A. Walmsley, Phys. Rev. A 86, 043820 (2012).
  * Hu  _et al._ (2016) L.-Y. Hu, J.-N. Wu, Z. Liao, and M. S. Zubairy, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 175504 (2016).
  * Hu  _et al._ (2017) L. Hu, Z. Liao, and M. S. Zubairy, Phys. Rev. A 95, 012310 (2017).
  * Glauber (1963) R. J. Glauber, Phys. Rev. 130, 2529 (1963).
  * Veitch  _et al._ (2014) V. Veitch, S. H. Mousavian, D. Gottesman, and J. Emerson, New J. Phys. 16, 013009 (2014).
  * Fannes  _et al._ (1992) M. Fannes, B. Nachtergaele, and R. F. Werner, Commun. Math. Phys. 144, 443 (1992).
  * Frisk Kockum  _et al._ (2019) A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys 1, 19 (2019).
  * Periwal  _et al._ (2021) A. Periwal, E. S. Cooper, P. Kunkel, J. F. Wienand, E. J. Davis, and M. Schleier-Smith, Nat. 600, 630 (2021).
  * de Buy Wenniger  _et al._ (2022) I. M. de Buy Wenniger, S. E. Thomas, M. Maffei, S. C. Wein, M. Pont, A. Harouri, A. Lemaître, I. Sagnes, N. Somaschi, A. Auffèves, and P. Senellart, arXiv:2202.01109 (2022).
  * Ballester  _et al._ (2012) D. Ballester, G. Romero, J. J. García-Ripoll, F. Deppe, and E. Solano, Phys. Rev. X 2, 021007 (2012).
  * Blais  _et al._ (2021) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * Leibfried  _et al._ (2003) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Rev. Mod. Phys. 75, 281 (2003).
  * Lv  _et al._ (2018) D. Lv, S. An, Z. Liu, J.-N. Zhang, J. S. Pedernales, L. Lamata, E. Solano, and K. Kim, Phys. Rev. X 8, 021027 (2018).
  * Gerry and Knight (2005) C. Gerry and P. Knight,  _Introductory quantum optics_ (Cambridge university press, 2005).
  * Haroche and Raimond (2006) S. Haroche and J.-M. Raimond,  _Exploring the quantum: atoms, cavities, and photons_ (Oxford university press, 2006).

##  Appendix A Jaynes-Cummings model

The first part of this appendix provides a brief overview of the eigenproblem
solution for the Jaynes-Cummings model (for further information, see
references Gerry and Knight (2005); Haroche and Raimond (2006)). In the
following, we derive the final reduced states of the cavity S and the atom C,
and explicitly determine the set of atomic states that satisfy the catalytic
constraint [Eq. (1)].

###  A.1 Eigenproblem

The Jaynes-Cummings interaction Hamiltonian
$H_{\text{int}}=g\left(a\sigma_{+}+a^{\dagger}\sigma_{-}\right)$ couples pairs
of atom-field states $\\{\ket{n+1,g},\ket{n,e}\\}$. Consequently, the
Hamiltonian $H$ decouples into a direct product of $2\times 2$-matrix
Hamiltonians, i.e., $H=\bigoplus_{n=0}^{\infty}H^{(n)}$, where

$H^{(n)}\begin{bmatrix}\ket{n+1,g}\\\
\ket{n,e}\end{bmatrix}=\begin{bmatrix}(n+1/2)\omega&g\sqrt{n+1}\,\\\ \\\
g\sqrt{n+1}&(n+1/2)\omega\end{bmatrix}\begin{bmatrix}\ket{n+1,g}\\\
\ket{n,e}\end{bmatrix}.$ (16)

The eigenvalue problem for this Hamiltonian yields the eigenfrequencies

$\omega^{(n)}_{\pm}=\quantity(n+\frac{1}{2})\omega\pm\frac{1}{2}\mu_{n},$ (17)

where $\mu_{n}=2g\sqrt{n+1}$ is the $n$-photon Rabi frequency. In the
resonance regime, the corresponding eigenstates are:

$\displaystyle\ket{n,+}=\frac{1}{\sqrt{2}}(\ket{n+1,g}+\ket{n,e}),\qquad\ket{n,-}=\frac{1}{\sqrt{2}}(\ket{n+1,g}-\ket{n,e}).$
(18)

The time evolution operator takes the form of

$\displaystyle U(t)=$ $\displaystyle e^{i\omega
t}\left|0,g\right\rangle\\!\\!\left\langle
0,g\right|+\sum_{n=0}^{\infty}e^{-i\quantity(n+\frac{1}{2})\omega
t}\quantity{\cos\frac{\mu_{n}t}{2}\Bigl{(}\left|n+1,g\right\rangle\\!\\!\left\langle
n+1,g\right|+\left|n,e\right\rangle\\!\\!\left\langle
n,e\right|\Bigl{)}-i\sin\frac{\mu_{n}t}{2}\Bigl{(}\left|n+1,g\right\rangle\\!\\!\left\langle
n,e\right|+\left|n,e\right\rangle\\!\\!\left\langle n+1,g\right|\Bigl{)}}.$
(19)

###  A.2 Reduced states of subsystems

The unitary operator given by equation (19) describes the dynamics of the
Jaynes-Cummings model in resonance and under the rotating-wave approximation.
After the interaction, the joint system
$\sigma_{\textsf{S}C}(t)=U(t)(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})U^{\dagger}$
becomes correlated. Let the cavity and the atom be prepared in general mixed
states, i.e.
$\rho_{\textsf{S}}=\sum^{\infty}_{n,m=0}p_{n,m}\left|n\right\rangle\\!\\!\left\langle
m\right|$ and $\chi_{\textsf{C}}=q\left|g\right\rangle\\!\\!\left\langle
g\right|+r\left|g\right\rangle\\!\\!\left\langle
e\right|+r^{*}\left|e\right\rangle\\!\\!\left\langle
g\right|+(1-q)\left|e\right\rangle\\!\\!\left\langle e\right|$. Then, the
reduced state of the cavity at time $t$ is obtained by taking the partial
trace over the atom’s degrees of freedom, i.e.
$\sigma_{\textsf{S}}(t):=\text{tr}_{\textsf{C}}[{\sigma_{\textsf{S}C}}(t)]$.
More specifically,

$\displaystyle\sigma_{\textsf{S}}(t)=$ $\displaystyle
qp_{0,0}\outerproduct{0}{0}+\sum_{n=0}^{\infty}e^{i(n+1)\omega
t}\quantity[qp_{0,n+1}\cos\quantity(\frac{\mu_{n}t}{2})+irp_{0,n}\sin\quantity(\frac{\mu_{n}t}{2})]\left|0\right\rangle\\!\\!\left\langle
n+1\right|+\text{h.c}$ $\displaystyle+\sum_{n,m=0}^{\infty}e^{-i(n-m)\omega
t}qp_{n+1,m+1}\quantity[\sin\quantity(\frac{\mu_{n}t}{2})\sin\quantity(\frac{\mu_{m}t}{2})\left|n\right\rangle\\!\\!\left\langle
m\right|+\cos\quantity(\frac{\mu_{n}t}{2})\cos\quantity(\frac{\mu_{m}t}{2})\left|n+1\right\rangle\\!\\!\left\langle
m+1\right|]$ $\displaystyle+\sum_{n,m=0}^{\infty}e^{-i(n-m)\omega
t}(1-q)p_{n,m}\quantity[\cos\quantity(\frac{\mu_{n}t}{2})\cos\quantity(\frac{\mu_{m}t}{2})\left|n\right\rangle\\!\\!\left\langle
m\right|+\sin\quantity(\frac{\mu_{n}t}{2})\sin\quantity(\frac{\mu_{m}t}{2})\left|n+1\right\rangle\\!\\!\left\langle
m+1\right|]$ $\displaystyle+\sum_{n,m=0}^{\infty}ie^{-i(n-m)\omega
t}\cos\quantity(\frac{\mu_{n}t}{2})\sin\quantity(\frac{\mu_{m}t}{2})\quantity[p_{n+1,m}r\outerproduct{n+1}{m+1}+p_{n,m+1}r^{*}\outerproduct{n}{m}]+\text{h.c}.$
(20)

The atomic state is obtained by marginalizing over the photonic degrees of
freedom, i.e. $\chi(t):=\text{tr}_{\textsf{S}}[{\sigma_{\textsf{S}C}}(t)]$,
which leads to

$\displaystyle\chi_{\textsf{C}}(t)=q(t)\left|g\right\rangle\\!\\!\left\langle
g\right|+r(t)\left|g\right\rangle\\!\\!\left\langle
e\right|+r^{*}(t)\left|e\right\rangle\\!\\!\left\langle
g\right|+[1-q(t)]\left|e\right\rangle\\!\\!\left\langle e\right|,$ (21)

Note that we identify $r:=r(0)$ and $q:=q(0)$ in what follows. The
coefficients $q(t)$ and $r(t)$ are given by

$\displaystyle q(t)$
$\displaystyle=q\sum_{n=0}^{\infty}p_{n}\cos^{2}\quantity(\frac{\mu_{n-1}t}{2})+(1-q)p_{n}\sin^{2}\quantity(\frac{\mu_{n}t}{2})+\textrm{Re}[irp_{n+1,n}]\sin{\quantity(\mu_{n}t)},$
(22) $\displaystyle r(t)$ $\displaystyle=-ie^{i\omega
t}\sum_{n=0}^{\infty}p_{n,n+1}\sin\quantity(\frac{\mu_{n}t}{2})\cos\quantity(\frac{\mu_{n+1}t}{2})+re^{i\omega
t}\sum_{n=0}^{\infty}p_{n}\cos\quantity(\frac{\mu_{n-1}t}{2})\cos\quantity(\frac{\mu_{n}t}{2})+r^{*}\sum_{n=0}^{\infty}p_{n,n+2}\sin\quantity(\frac{\mu_{n}t}{2})\sin\quantity(\frac{\mu_{n+1}}{2})$
$\displaystyle\hskip 224.7766pt+ie^{i\omega
t}q\sum_{n=0}^{\infty}p_{n,n+1}\sin\quantity(\frac{\mu_{n}t}{2})\quantity[\cos\quantity(\frac{\mu_{n-1}t}{2})+\cos\quantity(\frac{\mu_{n+1}t}{2})].$
(23)

###  A.3 Catalytic constraint

Here we determine the set of atomic states that evolve catalytically by
explicitly solving the catalytic constraint from Eq. (1). More specifically,
we are looking for the solution to the following operator equation:

$\displaystyle\chi_{\textsf{C}}(\tau)=\Tr_{\textsf{S}}\\{U(\tau)[\rho_{\textsf{S}}\otimes\chi_{\textsf{C}}(\tau)]U(\tau)^{\dagger}\\},$
(24)

for a fixed time $\tau$. From this point forward, we will abbreviate the
diagonal elements of the state $\rho_{\textsf{S}}$ as $p_{n}:=p_{n,n}$. To
obtain the set of states that satisfy Eq.(24), we first define the auxiliary
functions:

$\displaystyle\begin{split}\tilde{a}_{1}(t)&=e^{i\omega
t}\sum_{n=0}^{\infty}p_{n}\cos(gt\sqrt{n})\cos(gt\sqrt{n+1})-1,\\\
\tilde{a}_{2}(t)&=ie^{i\omega
t}\sum_{n=0}^{\infty}p_{n,n+1}\sin(gt\sqrt{n+1})\quantity[\cos(gt\sqrt{n})+\cos(gt\sqrt{n+2})],\\\
\tilde{a}_{3}(t)&=e^{i\omega
t}\sum_{n=0}^{\infty}p_{n,n+2}\sin(gt\sqrt{n+1})\sin(gt\sqrt{n+2}),\\\
\tilde{a}_{4}(t)&=-ie^{i\omega
t}\sum_{n=0}^{\infty}p_{n,n+1}\sin(gt\sqrt{n+1})\cos(gt\sqrt{n+2}).\end{split}$
(25)

Note that $\tilde{a}_{i}=e^{i\omega t}a_{i}$. Next, we observe that Eq. (24)
gives rise to a set of two equations with two variables. By considering the
ground state occupation $q(t)$, we find that the states satisfying Eq. (24)
are given by:

$q(t)=\frac{\sum_{n=0}^{\infty}p_{n}\sin^{2}(gt\sqrt{n+1})+\operatorname{Re}\quantity[ir(t)p_{n+1,n}]\sin(2gt\sqrt{n+1})}{\sum_{n=0}^{\infty}(p_{n}+p_{n+1})\sin^{2}(gt\sqrt{n+1})},$
(26)

whereas the coherence $r(t)$ obeys the equation

$r(t)\tilde{a}_{1}(t)+q(t)\tilde{a}_{2}(t)+r^{*}(t)\tilde{a}_{3}(t)+\tilde{a}_{4}(t)=0,$
(27)

whose solution is given by

$r(t)=\frac{\tilde{a}_{3}(t)\tilde{a}^{*}_{4}(t)-\tilde{a}^{*}_{1}(t)\tilde{a}_{4}(t)}{|\tilde{a}_{1}(t)|^{2}-|\tilde{a}_{3}(t)|^{2}}+\frac{\tilde{a}_{3}(t)\tilde{a}^{*}_{2}(t)-\tilde{a}^{*}_{1}(t)\tilde{a}_{2}(t)}{|\tilde{a}_{1}(t)|^{2}-|\tilde{a}_{3}(t)|^{2}}q(t).$
(28)

Substituting Eq. (28) into Eq. (26), we find that

$q(t)=\frac{\sum_{n=0}^{\infty}p_{n}\sin^{2}(gt\sqrt{n+1})+\operatorname{Re}\quantity[i\quantity(\frac{\tilde{a}_{3}(t)\tilde{a}^{*}_{4}(t)-\tilde{a}^{*}_{1}(t)\tilde{a}_{4}(t)}{|\tilde{a}_{1}(t)|^{2}-|\tilde{a}_{3}(t)|^{2}})p_{n+1,n}]\sin(2gt\sqrt{n+1})}{\sum_{n=0}^{\infty}(p_{n+1}+p_{n})\sin^{2}(gt\sqrt{n+1})]-\operatorname{Re}\quantity[i\quantity(\frac{\tilde{a}_{3}(t)\tilde{a}^{*}_{2}(t)-\tilde{a}^{*}_{1}(t)\tilde{a}_{2}(t)}{|\tilde{a}_{1}(t)|^{2}-|\tilde{a}_{3}(t)|^{2}})p_{n+1,n}]\sin(2gt\sqrt{n+1})}.$
(29)

Therefore, for a given value of $g$ and time $\tau$, Eqs. (28) and (29)
uniquely determine a state of the catalyst.

##  Appendix B Quantum catalysis and higher moments of observables

In this Appendix we extend the analysis of a catalytic evolution presented in
Sec. III. In particular, we derive an expression for the $k$-th moment of an
arbitrary observable $O_{\textsf{S}}$ on the system S undergoing a catalytic
evolution. Specifying this expression to the second moment $(k=2)$ and the
particle number observable $(O_{\textsf{S}}=n_{\textsf{S}})$ leads to Eq. (9)
as stated in the main text. Then, using the results of Appendix A, we obtain
an explicit expression for the second moment in a catalytic evolution as
specified by the Jaynes-Cummings Hamiltonian.

###  B.1 Higher moments of observables under catalytic evolution

Let
$O_{\textsf{SC}}=O_{\textsf{S}}\otimes\mathbbm{1}_{\textsf{C}}+\mathbbm{1}_{\textsf{S}}+O_{\textsf{C}}$
be and additive and conserved observable on the joint system SC. In a
catalytic protocol governed by a unitary $U$ we have
$\sigma_{\textsf{SC}}=U(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})U^{\dagger}$
such that $\sigma_{\textsf{C}}=\chi_{\textsf{C}}$. The conservation assumption
means that $[U,O_{\textsf{SC}}]=0$. The $k$-th moment of observable $A$ on the
system S before and after the catalytic process are respectively given by

$\displaystyle\langle
O^{k}_{\textsf{S}}\rangle_{\rho}:=\Tr[O_{\textsf{S}}^{k}\rho_{\textsf{S}}],\qquad\langle
O^{k}_{\textsf{S}}\rangle_{\sigma}:=\Tr[O_{\textsf{S}}^{k}\sigma_{\textsf{S}}]=\Tr[(O_{\textsf{S}}^{k}\otimes\mathbbm{1}_{\textsf{C}})\sigma_{\textsf{SC}}].$
(30)

Let us write
$O_{\textsf{SC}}^{k}=O_{\textsf{S}}^{k}\otimes\mathbbm{1}_{\textsf{C}}+\mathbbm{1}_{\textsf{C}}\otimes
O_{\textsf{C}}^{k}+\Delta_{k}$, where we defined an auxiliary observable

$\displaystyle\Delta_{k}:=O_{\textsf{SC}}^{k}-O_{\textsf{S}}^{k}\otimes\mathbbm{1}_{\textsf{C}}-\mathbbm{1}_{\textsf{C}}\otimes
O_{\textsf{C}}^{k}=\sum_{i=1}^{k-1}\binom{k-1}{i}O^{k-i-1}_{\textsf{S}}\otimes
O^{i}_{\textsf{C}}.$ (31)

For example, for $k=2$ we have $\Delta_{2}=2O_{\textsf{S}}\otimes
O_{\textsf{C}}$. Using this to express $\langle
O^{k}_{\textsf{S}}\rangle_{\sigma}$ yields

$\displaystyle\langle
O^{k}_{\textsf{S}}\rangle_{\sigma}=\Tr[O_{\textsf{SC}}^{k}\sigma_{\textsf{SC}}]-\Tr[(\mathbbm{1}_{\textsf{S}}\otimes
O_{\textsf{C}}^{k})\sigma_{\textsf{SC}}]-\Tr[\Delta_{k}\sigma_{\textsf{SC}}]=\Tr[O_{\textsf{SC}}^{k}\sigma_{\textsf{SC}}]-\langle
O_{\textsf{C}}^{k}\rangle_{\sigma}-\Tr[\Delta_{k}\sigma_{\textsf{SC}}].$ (32)

Importantly, due to the conservation law $[U,O_{\textsf{SC}}]=0$, we have
$\Tr[O_{\textsf{SC}}^{k}\sigma_{\textsf{SC}}]=\Tr[O_{\textsf{SC}}^{k}\rho_{\textsf{SC}}]=\langle
O_{\textsf{S}}^{k}\rangle_{\rho}+\langle
O_{\textsf{C}}^{k}\rangle_{\rho}+\Tr[\Delta_{k}{\rho_{\textsf{SC}}}]$.
Therefore we can rewrite Eq. (32) as

$\displaystyle\langle O^{k}_{\textsf{S}}\rangle_{\sigma}=\langle
O_{\textsf{S}}^{k}\rangle_{\rho}+\langle
O_{\textsf{C}}^{k}\rangle_{\rho}-\langle
O_{\textsf{C}}^{k}\rangle_{\sigma}+\Tr[\Delta_{k}(\rho_{\textsf{SC}}-\sigma_{\textsf{SC}})]=\langle
O_{\textsf{S}}^{k}\rangle_{\rho}+\Tr[\Delta_{k}(\rho_{\textsf{SC}}-\sigma_{\textsf{SC}})],$
(33)

where we used the fact that $\langle O_{\textsf{C}}^{k}\rangle_{\rho}=\langle
O_{\textsf{C}}^{k}\rangle_{\sigma}$ due to the catalytic constraint. By taking
$k=2$ we obtain the desired result.

###  B.2 Second moment of photon statistics in the catalytic Jaynes-Cummings
evolution

Consider the second moment of photon statistics

$\displaystyle\langle n^{2}_{\textsf{S}}\rangle_{\sigma}=\langle
n^{2}_{\textsf{S}}\rangle_{\rho}+2\left[(1-q)\,\langle
n_{\textsf{S}}\rangle_{\rho}-\langle
n_{\textsf{S}}\otimes\outerproduct{e}{e}_{\textsf{C}}\rangle_{\sigma}\right],$
(34)

where $q:=\langle e|\chi_{\textsf{C}}|e\rangle$ is the excited-state
occupation of the catalyst,
$\sigma=U(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})U^{\dagger}$ and
$\Tr_{\textsf{S}}[\sigma]=\chi_{\textsf{C}}$. Let us focus on the following
term:

$\displaystyle\langle
n_{\textsf{S}}\otimes\outerproduct{e}{e}_{\textsf{C}}\rangle_{\sigma}$
$\displaystyle=\Tr\left[U^{\dagger}(n_{\textsf{S}}\otimes\outerproduct{e}{e}_{\textsf{C}})U(\rho_{\textsf{S}}\otimes\chi_{\textsf{C}})\right]=\sum_{k=0}^{\infty}k\Tr[U^{\dagger}\outerproduct{k,e}{k,e}U(\rho_{\textsf{S}}\otimes\chi)],$
(35)

Using Eq. (19), we can write

$\displaystyle U^{\dagger}\ket{k,e}$ $\displaystyle=e^{i(k+\frac{1}{2})\chi
t}\left(c_{k}\ket{k,e}-is_{k}\ket{k+1,g}\right),$ (36)

where $c_{k}:=\cos(gt\sqrt{k+1})$ and $s_{k}:=\sin(gt\sqrt{k+1})$.
Substituting the above result into Eq. (35) leads to

$\displaystyle\langle
n_{\textsf{S}}\otimes\outerproduct{e}{e}_{\textsf{C}}\rangle_{\sigma}=\sum_{k=0}^{\infty}k\left[(1-q)c_{k}^{2}p_{k,k}+2s_{k}c_{k}\text{Im}(p_{k+1,k}r)+(1-q)s_{k}^{2}p_{k+1,k+1}\right].$
(37)

##  Appendix C Exploring catalytic evolution in the Jaynes-Cummings model

In this Appendix we examine different aspects of the catalytic evolution in
the Jaynes-Cummings model.

###  C.1 Which states of the cavity cannot have their non-classicality
enhanced?

Let us start by examining Eqs. (28) and (29) under the assumption that the
initial state of the cavity is an incoherent mixture of Fock states, i.e.
$\rho_{\textsf{S}}=\sum_{n,m}p_{n,m}\outerproduct{n}{m}$ with $p_{n,m}=0$ if
$n\neq m$. In this case, from Eqs. (13) and (14) we can infer that the only
feasible states of the catalyst are those with $q_{\text{coh}}=0$ and $r=0$.
Consequently, the atomic state is incoherent in the energy basis and its
ground state occupation takes the form of

$q=\frac{\sum_{n=0}^{\infty}p_{n}\sin^{2}(gt\sqrt{n+1})}{\sum_{n=0}^{\infty}(p_{n}+p_{n+1})\sin^{2}(gt\sqrt{n+1})}.$
(38)

Let us now demonstrate that the second-order auto-correlation function
$g^{(2)}$ of a pure Fock state, i.e.
$\rho_{\textsf{S}}=\outerproduct{k}{k}_{\textsf{S}}$, cannot decrease under a
catalytic evolution. To do this, we will prove that the quantity

$g^{(2)}(\sigma_{\textsf{S}})-g^{(2)}(\rho_{\textsf{S}})=\frac{\Tr[(a^{\dagger
2}a^{2})\sigma_{\textsf{S}}]}{\Tr[(a^{\dagger}a)\sigma_{\textsf{S}}]}-\frac{\Tr[(a^{\dagger
2}a^{2})\rho_{\textsf{S}}]}{\Tr[(a^{\dagger}a)\rho_{\textsf{S}}]}=\frac{\Tr[(a^{\dagger
2}a^{2})(\sigma_{\textsf{S}}-\rho_{\textsf{S}})]}{\Tr[(a^{\dagger}a)\rho_{\textsf{S}}]}$
(39)

is nonnegative. This will be accomplished by showing that the following
inequality holds:

$\Tr[(a^{\dagger 2}a^{2})(\sigma_{\textsf{S}}-\rho_{\textsf{S}})]\geq 0.$ (40)

In what follows, we will omit the index S as well as any explicit reference to
the variables’ dependence on time. When the initial state of the cavity is a
Fock state, then Eq. (A.2) implies that the state of the cavity after the
catalytic evolution is given by

$\displaystyle\\!\\!\\!\sigma=\quantity[(1-q)\cos^{2}\quantity(gt\sqrt{k+1})+q\cos^{2}\quantity(gt\sqrt{k})]\outerproduct{k}{k}+q\sin^{2}\quantity(gt\sqrt{k})\outerproduct{k-1}{k-1}+(1-q)\sin^{2}\quantity(gt\sqrt{k+1})\outerproduct{k+1}{k+1},$
(41)

where $q$ is determined by Eq. (38), which for this particular case takes the
form:

$q=\frac{\sin^{2}\quantity(gt\sqrt{k+1})}{\sin^{2}\quantity(gt\sqrt{k+1})+\sin^{2}\quantity(gt\sqrt{k})}\quad\text{and}\quad
1-q=\frac{\sin^{2}\quantity(gt\sqrt{k})}{\sin^{2}\quantity(gt\sqrt{k+1})+\sin^{2}\quantity(gt\sqrt{k})}.$
(42)

By substituting Eq.(42) into Eq.(A.2) and introducing the notation
$\psi=\frac{1}{2}(\outerproduct{n-1}{n-1}+\outerproduct{n+1}{n+1})$, we
obtain:

$\sigma=\frac{\sin^{2}\quantity(gt\sqrt{k})\cos^{2}\quantity(gt\sqrt{k+1})+\sin^{2}\quantity(gt\sqrt{k+1})\cos^{2}\quantity(gt\sqrt{k})}{\sin^{2}\quantity(gt\sqrt{k+1})+\sin^{2}\quantity(gt\sqrt{k})}\rho+\frac{2\sin^{2}\quantity(gt\sqrt{k})\sin^{2}\quantity(gt\sqrt{k+1})}{\sin^{2}\quantity(gt\sqrt{k+1})+\sin^{2}\quantity(gt\sqrt{k})}\psi.$
(43)

Alternatively, we can express the above equation as $\sigma=t\rho+(1-t)\psi$,
where

$t=\frac{\sin^{2}\quantity(gt\sqrt{k})\cos^{2}\quantity(gt\sqrt{k+1})+\sin^{2}\quantity(gt\sqrt{k+1})\cos^{2}\quantity(gt\sqrt{k})}{\sin^{2}\quantity(gt\sqrt{k+1})+\sin^{2}\quantity(gt\sqrt{k})}.$
(44)

With these results at hand, we can manipulate Eq. (40) to obtain

$\displaystyle\Tr[(a^{\dagger 2}a^{2})(\sigma-\rho)]$
$\displaystyle=(1-t)\Tr[(a^{\dagger
2}a^{2})(\psi-\rho)]=\frac{(1-t)}{2}\Tr[(a^{\dagger
2}a^{2})(\outerproduct{k+1}{k+1}+\outerproduct{k-1}{k-1})]-(1-t)\Tr[(a^{\dagger
2}a^{2})\outerproduct{k}{k}]$
$\displaystyle=\frac{(1-t)}{2}\quantity[(k+1)k+(k-1)(k-2)-2k(k-1)]=(1-t)\geq
0.$ (45)

Thus, we conclude that second-order coherence $g^{(2)}$ can only increase
during a catalytic processes involving a pure Fock state, i.e.

$g^{(2)}(\sigma)-g^{(2)}(\rho)\geq 0.$ (46)

It is not clear whether there exist other states of the cavity whose non-
classicality, as witnessed by the $g^{(2)}$ function, cannot be further
enahnced. We leave this interesting problem for future work.

###  C.2 Which states of the atom allow to generate non-classicality?

An interesting question is which states of the atom can satisfy the catalytic
constraint from Eq. (1) for some value of the stopping time $\tau$, and for a
fixed initial (coherent) state of the cavity, i.e.
$\rho_{\textsf{S}}=\outerproduct{\alpha}{\alpha}$ with amplitude
$|\alpha|^{2}$. In what follows we will refer to the set of all atomic states
that satisfy Eq. (1) as the _set of catalytic states_. Moreover, we observe
that not every state in this set can lead to a final state of the cavity that
is non-classical.

To investigate these questions, we characterize numerically the catalytic set
for three initial coherent states of the cavity corresponding to
$\alpha\in\\{0.2,1/\sqrt{2},25\\}$. Results are presented in Fig. 5, where the
catalytic set is shown in grey. Additionally, we also determine which states
in the catalytic set can generate non-classicality, focusing on the $g^{(2)}$
function here. Similarly to Fig. 4(b) of the main text, we represent these
states in colour, the latter indicating the level of non-classicality being
generated. Interestingly, these sets vary significantly with $|\alpha|$. Also,
when $|\alpha|$ is large, catalytic states are distributed close to the
equatorial plane of the Bloch sphere, and we could find no instance where non-
classicality is generated.

Figure 5: Set of catalytic state for different initial state preparation.
Atomic states (projection of the $y-z$ plane of the Bloch ball) that satisfy
the catalytic constraint (gray) for (a) $\alpha=0.2$, (b) $\alpha=1/\sqrt{2}$,
and (c) $\alpha=25$ are highlighted by a uniform gradient if they produce
$g^{2}(\sigma_{\textsf{S}})<1$ All panels were generated by imposing a limit
on interaction time with $\tau\leq 100$ and taking $10^{6}$ samples. The
parameters are $\omega=2\pi$ and $g=\pi$.

###  C.3 Boosting non-classicality via a catalytic process for incoherent
mixtures of Fock states

For an incoherent mixture of Fock states
$\rho=\sum_{n}p_{n}\outerproduct{n}{n}$, the second-order coherence is given
by

$g^{(2)}(\rho)=\frac{\sum_{n=0}^{\infty}n(n-1)p_{n}}{(\sum_{n=0}^{\infty}np_{n})^{2}}.$
(47)

Assuming that the initial state of the cavity is prepared in a state
$\rho_{\textsf{S}}=\frac{1}{4}\outerproduct{0}{0}+\frac{3}{4}\outerproduct{2}{2})$,
then its second-order coherence is $g^{(2)}(\rho_{\textsf{S}})=2/3$. According
to Eq. (A.2), the state of the cavity after the catalytic protocol (at time
$t=\tau$) takes the form of

$\displaystyle\sigma_{\textsf{S}}=\frac{1}{4}\quantity[q+(1-q)\cos^{2}gt]\outerproduct{0}{0}+\frac{3q}{4}\quantity[\sin^{2}\quantity(g\tau\sqrt{2})+\frac{(1-q)}{4}\sin^{2}gt]\outerproduct{1}{1}$
$\displaystyle+\frac{3}{4}\quantity[q\cos^{2}\quantity(g\tau\sqrt{2})+(1-q)\cos^{2}\quantity(g\tau\sqrt{3})]\outerproduct{2}{2}$
$\displaystyle+\frac{3(1-q)}{4}\sin^{2}\quantity(g\tau\sqrt{3})\outerproduct{3}{3}$
(48)

where $q$ is determined by Eq. (38):

$q=\frac{\sin^{2}\quantity(g\tau)+3\sin^{2}\quantity(g\tau\sqrt{3})}{\sin^{2}\quantity(g\tau)+3\quantity[\sin^{2}\quantity(g\tau\sqrt{3})+\sin^{2}\quantity(g\tau\sqrt{2})]}.$
(49)

Using Eq. (47), the second-order coherence for the final state is given by

$g^{(2)}(\sigma_{\textsf{S}})=\frac{2}{3}\quantity{q\cos^{2}\quantity(g\tau\sqrt{2})+(1-q)\quantity[1+2\sin^{2}\quantity(g\tau\sqrt{3})]}.$
(50)

Therefore, for $g\tau=7.5\pi$, we obtain $g^{(2)}\approx 0.505$, indicating
that non-classicality in the mode was catalytically increased.

Interestingly, when the cavity is initially prepared in a mixture of Fock
states, we oberve that the catalyst must necessarily be in an incoherent state
for satisfying the catalytic constraint.


**footnotetext: Equal contribution

#  Improving Trustworthiness of AI Disease Severity Rating in Medical Imaging
with Ordinal Conformal Prediction Sets

Charles Lu*, Anastasios N. Angelopoulos*, Stuart Pomerantz

###### Abstract

The regulatory approval and broad clinical deployment of medical AI have been
hampered by the perception that deep learning models fail in unpredictable and
possibly catastrophic ways. A lack of statistically rigorous uncertainty
quantification is a significant factor undermining trust in AI results. Recent
developments in distribution-free uncertainty quantification present practical
solutions for these issues by providing reliability guarantees for black-box
models on arbitrary data distributions as formally valid finite-sample
prediction intervals. Our work applies these new uncertainty quantification
methods — specifically conformal prediction — to a deep-learning model for
grading the severity of spinal stenosis in lumbar spine MRI. We demonstrate a
technique for forming ordinal prediction sets that are guaranteed to contain
the correct stenosis severity within a user-defined probability (confidence
interval). On a dataset of 409 MRI exams processed by the deep-learning model,
the conformal method provides tight coverage with small prediction set sizes.
Furthermore, we explore the potential clinical applicability of flagging cases
with high uncertainty predictions (large prediction sets) by quantifying an
increase in the prevalence of significant imaging abnormalities (e.g. motion
artifacts, metallic artifacts, and tumors) that could degrade confidence in
predictive performance when compared to a random sample of cases.

##  1 Introduction

Although many studies have demonstrated high overall accuracy in automating
medical imaging diagnosis with deep-learning AI models, translation to actual
clinical deployment has proved difficult. It is widely observed that deep
learning algorithms can fail in bizarre ways and with misplaced confidence
[23, 14]. A core problem is a lack of _trust_ — a survey of radiologists found
they that although they thought AI tools add value to their clinical practice,
they would not trust AI for autonomous clinical use due to perceived and
experienced unreliability [1].

Herein, we present methods for endowing arbitrary AI systems with formal
mathematical guarantees that give clinicians explicit assurances about an
algorithm’s overall performance and most importantly for a given study. These
guarantees are _distribution-free_ —they work for any (pre-trained) model, any
(possibly unknown) data distribution, and in finite samples. Although such
guarantees do not solve the issue of trust entirely, the precise and formal
understanding of their model’s predictive uncertainty enables a clinician to
potentially work more assuredly in concert with AI assistance.

We demonstrate the utility of our methods using an AI system developed to
assist radiologists in the grading of spinal stenosis in lumbar MRI.
Degenerative spinal stenosis is the abnormal narrowing of the spinal canal
that compresses the spinal cord or nerve roots, often resulting in debility
from pain, limb weakness, and other neurological disorders. It is a highly
prevalent condition that affects working-age and elderly populations and
constitutes a heavy societal burden not only in the form of costly medical
care but from decreased workplace productivity, disability, and lowered
quality of life. The formal interpretation of spinal stenosis imaging remains
a challenging and time-consuming task even for experienced subspecialty
radiologists due to the complexity of spinal anatomy, pathology, and the MR
imaging modality. Using a highly accurate AI model to help assess the severity
of spinal stenosis on MRI could lower interpretation time and improve the
consistency of grading. [12] Yet, given the practical challenges of medical
imaging that can degrade model performance in any given exam, the adoption of
such tools will be low if clinicians encounter poor quality predictions
without a sense of when the model is more or less reliable. To bolster such
trust, we apply conformal prediction to algorithmic disease severity
classification sets in order to identify higher uncertainty predictions that
might merit special attention by the radiologist. Our main contributions are
the following:

  1. 1.

We develop new distribution-free uncertainty quantification methods for
ordinal labels.

  2. 2.

To our knowledge, we are the first to apply distribution-free uncertainty
quantification to the results of an AI model for automated stenosis grading of
lumbar spinal MRI.

  3. 3.

We identify a correlation between high prediction uncertainty in individual
cases and the presence of potentially contributory imaging features as
evaluated by a neuroradiologist such as tumors, orthopedic hardware artifacts,
and motion artifacts.

Figure 1: The output of our conformal procedure on a lumbar spine MRI.

##  2 Methods

To formally describe the problem, let the input $X_{\rm test}\in\mathcal{X}$,
$\mathcal{X}=\mathbb{R}^{H\times W\times D}$ be an MR image and the ground
truth label $Y_{\rm test}\in\mathcal{Y}$, $\mathcal{Y}=\\{0,...,K-1\\}$ be an
ordinal value representing the severity of the disease (higher values
indicating greater severity). We are given a pre-trained model, $\hat{f}$,
that takes in images and outputs a probability distribution over severities;
for example, $\hat{f}$ may be a 3D convolutional neural network with a softmax
function. Assume we also have a calibration dataset,
$\big{\\{}(X_{i},Y_{i})\big{\\}}_{i=1}^{n}$, of data points that the model has
not seen before. This calibration data should be composed of pairs of MR
images and their associated disease severity labels drawn i.i.d. Given a new
MR image $X_{\rm test}$, the task is to predict the (unobserved) severity,
$Y_{\rm test}$. In the usual multi-class classification setting, the output is
the label with the highest estimated probability,
$\hat{Y}(x)=\underset{y\in\\{1,...,K\\}}{\arg\max}\hat{f}(x)_{y}$. However,
$\hat{Y}(X_{\rm test})$ may be wrong, either because the learned model
$\hat{f}$ does not learn the relationship between MR images and severities
properly or because there is intrinsic randomness in this relationship that
cannot be accounted for by any algorithm (i.e. aleatoric uncertainty).

Our goal is to rigorously quantify this uncertainty by outputting a set of
probable disease severities that is guaranteed to contain the ground truth
severity on average. These _prediction sets_ will provide _distribution-free_
probabilistic guarantees, _i.e._ ones that do not depend on the model or
distribution of the data.

#### Ordinal Adaptive Prediction Sets (Ordinal APS)

Our approach uses conformal prediction with a novel score function designed
for ordinal labels. In particular, each prediction set will always be a
contiguous set of severities, and for any user-specified error rate $\alpha$,
prediction sets will contain the true label with probability $1-\alpha$. The
reader can refer to [25] and [27] for similar algorithms and exposition.

#### An Oracle Method

Imagine we had oracle access to the true probability distribution over
severities $P\big{(}Y_{\rm test}\mid X_{\rm test}\big{)}$ with associated
density function $f(x)_{y}$. A reasonable goal might then be to pick the set
with the smallest size while still achieving _conditional coverage_.

######  Definition 1 (conditional coverage).

A predicted set of severities $\mathcal{T}(X_{\rm test})$ has conditional
coverage if it contains the true severity with $1-\alpha$ probability no
matter what MR image is input, i.e.,

$\small\mathbb{P}\left(Y_{\rm test}\in\mathcal{T}(x)\mid X_{\rm
test}=x\right)\geq 1-\alpha\text{, for all }x\in\mathcal{X}.$ (1)

The clinical benefit of conditional coverage is that it essentially achieves a
per-patient guarantee as opposed to one that is averaged across patients.
Under conditional coverage, the uncertainty sets will work equally well for
all possible subgroups such as subpopulations from different patient
demographics.

Ignoring tie-breaking, we can succinctly describe the oracle prediction set as
follows:

$\displaystyle\mathcal{T}^{(\rm optimal)}(x)$
$\displaystyle=\left[l^{*}(x),\;u^{*}(x)\right],\text{ where }$ (2)
$\displaystyle\big{(}l^{*}(x),\;u^{*}(x)\big{)}$
$\displaystyle=\underset{\begin{subarray}{c}(l,u)\in\mathcal{Y}^{2}\\\ l\leq
u\end{subarray}}{\arg\min}\Bigg{\\{}u-l:\sum\limits_{j=l}^{u}f(x)_{j}\geq
1-\alpha\Bigg{\\}}.$

This set, $\mathcal{T}^{(\rm optimal)}$, is the smallest that satisfies (1).
Ideally, we would compute $\mathcal{T}^{(\rm optimal)}$ exactly, but we do not
have access to $f$, only its estimator $\hat{f}$, which may be arbitrarily
bad.

#### Ordinal Adaptive Prediction Sets

Naturally, the next step is to plug in our estimate of the probabilities,
$\hat{f}$, to (2). However, because $\hat{f}$ may be wrong, we must calibrate
the resulting set with conformal prediction, yielding a marginal coverage
guarantee. Our procedure is illustrated graphically in the right plot of
Figure 2; it corresponds to greedily growing the set outwards from the
maximally likely predicted severity (i.e. “mild” stenosis in this example).

######  Definition 2 (marginal coverage).

A predicted set of severities $\mathcal{T}$ has marginal coverage if it
contains the true severity on average over new MRIs, i.e.,

$\small\mathbb{P}\left(Y_{\rm test}\in\mathcal{T}(X_{\rm test})\right)\geq
1-\alpha.$ (3)

Figure 2: Comparison of an example prediction set chosen by naive least-
ambiguous set-valued (LAC) classifiers [26] and Ordinal APS methods; notice
that LAC does not respect ordinality, which may result in an incongruous
prediction set.

Marginal coverage is weaker than conditional coverage since it holds only on
average over the entire population, so coverage may be worse or better for
some subgroups. While conditional coverage is, in general, impossible [7], we
can hope to approximate conditional coverage by defining a set similar to
$\mathcal{T}^{(\rm optimal)}$ that uses $\hat{f}$. To that end, define a
sequence of sets indexed by a threshold $\lambda$,

$\displaystyle\mathcal{T}_{\lambda}(x)$
$\displaystyle=\;[\hat{l}^{*}(x),\;\hat{u}^{*}(x)],\text{ where }$ (4)
$\displaystyle\big{(}\hat{l}^{*}(x),\;\hat{u}^{*}(x)\big{)}$
$\displaystyle=\underset{\begin{subarray}{c}(l,\;u)\;\in\;\mathcal{Y}^{2}\\\
l\;\leq\;u\end{subarray}}{\arg\min}\Bigg{\\{}u-l:\sum\limits_{j=l}^{u}\hat{f}(j|x)\geq\lambda\Bigg{\\}}.$

Notice that as $\lambda$ grows, the sets grow, meaning they are _nested_ in
$\lambda$:

$\lambda_{1}\leq\lambda_{2}\implies\forall
x,\;\;\mathcal{T}_{\lambda_{1}}(x)\subseteq\mathcal{T}_{\lambda_{2}}(x).$ (5)

The key is to pick a value of $\lambda$ such that the resulting set satisfies
(3). The following algorithm takes as input $\mathcal{T}_{\lambda}$ and
outputs our choice of $\lambda$:

$\small\mathcal{A}(\mathcal{T}_{\lambda};\;\alpha)=\inf\Bigg{\\{}\lambda:\sum\limits_{i=1}^{n}\mathbbm{1}\left\\{Y_{i}\in\mathcal{T}_{\lambda}(X_{i})\right\\}\;\geq\;\lceil(n+1)(1-\alpha)\rceil\Bigg{\\}}.$
(6)

The key to our guarantee is the quantity $\lceil(n+1)(1-\alpha)\rceil$, which
is slightly larger than the naive choice $n(1-\alpha)$ and helps us correct
for the model’s deficiencies; see [2] for details on this statistical
argument. Using this algorithm, approximated in Algorithm 1, results in a
marginal coverage guarantee.

######  Theorem 1 (Conformal coverage guarantee).

Let $(X_{1},Y_{1})$, $(X_{2},Y_{2})$, …, $(X_{n},Y_{n})$ and $(X_{\rm
test},Y_{\rm test})$ be an i.i.d. sequence of MRIs and paired severities and
let $\hat{\lambda}=\mathcal{A}(\mathcal{T}_{\lambda},\;\alpha)$. Then
$\mathcal{T}_{\hat{\lambda}}$ satisfies (3), i.e., it contains the true label
with probability $1-\alpha$.

This theorem holds for any data distribution or machine learning model, any
number of calibration data points $n$, and any possible sequence of nested
sets that includes $\mathcal{Y}$ (see the formal version and proof in Appendix
A).

#### Implementing Ordinal Adaptive Prediction Sets

In practice, Ordinal APS has two undesirable properties: computing
$\mathcal{T}_{\lambda}(x)$ requires a combinatorial search of the set of
possible severities, and $\mathcal{T}_{\lambda}(x)$ may not include the point
prediction $\hat{Y}$. In practice, we therefore approximate Ordinal APS
greedily as described below in Algorithm 1.

Algorithm 1 Pseudocode for approximately computing $\mathcal{T}_{\lambda}(x)$

Input: Parameter $\lambda$; underlying predictor $\hat{f}$; input
$x\in\mathcal{X}$.

Output: $\mathcal{T}_{\lambda}(x)$.

1:$\mathcal{T}_{\lambda}(x)\leftarrow\arg\max\hat{f}(x)$

2:$q\leftarrow 0$

3:while $q\leq\lambda$ do

4:
$S\leftarrow\\{\min\mathcal{T}_{\lambda}(x)-1,\;\max\mathcal{T}_{\lambda}(x)+1\\}$

5: $y\leftarrow\underset{y^{\prime}\in
S}{\arg\max}\,\hat{f}(x)\mathbbm{1}\left\\{y^{\prime}\in\\{1,...,K\\}\right\\}$

6: $q\leftarrow q+\hat{f}(x)_{y}$

7: $\mathcal{T}_{\lambda}(x)\leftarrow\mathcal{T}_{\lambda}(x)\;\cup\;\\{y\\}$

The algorithm always contains $\hat{Y}$ and requires only $\mathcal{O}(n)$
computations; furthermore, it usually results in exactly the same sets as the
exact method in our experiments, which have a small value of $K$.

Note that the approximate choice of $\mathcal{T}_{\lambda}(x)$ described in
Algorithm 1 is still nested, and thus we can still guarantee coverage (see
Corollary A.1 for a formal statement and proof).

##  3 Experiments

We compare Ordinal Adaptive Prediction Sets to two other conformal methods:
Least Ambiguous set-valued Classifier (LAC) [26] and Ordinal Cumulative
Distribution Function (CDF). LAC uses the softmax score of the true class as
the conformal score function. LAC theoretically gives the smallest average set
size but sacrifices conditional coverage to achieve this. Additionally, LAC
does not respect ordinality and thus may output non-contiguous prediction
sets, which are inappropriate in an ordinal setting such as disease severity
rating. The Ordinal CDF method starts at the highest prediction score and then
inflates the intervals by $\lambda$ in quantile-space; in that sense, it is
similar to a discrete version of conformalized quantile regression [15, 24].
We only use the non-randomized versions of these methods.

We evaluate these three conformal methods on a deep learning system previously
developed for automated lumbar spine stenosis grading in MRI, DeepSPINE [21].
The deep learning system consists of two convolutional neural networks – one
to segment out and label each vertebral body and disc-interspace and the other
to perform multi-class ordinal stenosis classification for three different
anatomical sites (the central canal and right and left neuroforamina) at each
intervertebral disc level for a total of up to 18 gradings per patient.

Figure 3: Example of a case without lumbar spinal stenosis and with severe
lumbar spinal stenosis.

For each MRI exam, the associated radiology report was automatically parsed
for keywords indicative of the presence and severity of stenosis for each of
the 6 vertebral disc levels (T12-L1 through L5-S1) to extract ground truth
labels for a total of 6,093 gradings. Each grading was assigned to a value on
a four-point ordinal scale of stenosis severity: 0 (no stenosis), 1 (mild
stenosis), 2 (moderate stenosis), and 3 (severe stenosis). Examples of
patients with and without stenosis are shown in Figure 3.

For our experiments, we treat the stenosis grading model as a static model and
only use it for inference to make predictions. We then process these predicted
scores with the split-conformal techniques described in Section 2. This
scenario would most closely reflect the clinical reality of incorporating
regulated, third-party AI software medical devices, which would likely not
permit users access to the AI model beyond the ability to make predictions.

Our code and analysis used in the quantitative experiments are made available
here: https://github.com/clu5/lumbar-conformal.

###  3.1 Quantitative Experiments

We use the predicted softmax scores from a held-out set of MRI exams from 409
patients, comprising 6,093 disc level stenosis severity predictions to
calibrate and evaluate each conformal methods. We randomly include 5% of
patients in the calibration set and reserve the remainder for evaluating
coverage and set size. We evaluate performance at several different $\alpha$
thresholds, $\alpha\in\\{0.2,0.15,0.1,0.05,0.01\\}$, and average results over
100 random trials.

Figure 4: Empirical coverage and set size for three conformal prediction
methods for $\alpha\in\\{0.2,0.15,0.1,0.05,0.01\\}$ (averaged over 100 trials
and shown with $\pm$ 1 standard deviation).

As expected, all three conformal methods empirically achieve the desired
marginal coverage as guaranteed by Theorem 1. However, Ordinal CDF requires a
much larger set size to attain proper coverage than either Naive LAC or
Ordinal APS for all values of $\alpha$ (shown in Figure 4).

In addition, while aggregate coverage is satisfied for each method, we find
significant differences in class-conditional coverage (i.e. prediction sets
stratified by the true stenosis severity label), which is shown in Figure 5.
We see that prediction sets for “mild stenosis” and “moderate stenosis” grades
have lower average coverage and larger set sizes than prediction sets for “no
stenosis” and “severe stenosis” grades. These differences may be partly
attributed to the fact that the “no stenosis” class constitutes the majority
of the label distribution (67%) and so may be easier for a model trained on
this distribution to classify. Additionally, “mild stenosis” and “moderate
stenosis” grades may be more challenging to differentiate than “no stenosis”
and “severe stenosis” grades, reflecting the greater inherent variability and
uncertainty in the ground truth ratings.

Figure 5: Comparing coverage and set size when stratified by ground-truth
stenosis severity grading. Figure 6: Comparing coverage and count of number of
sets with a particular size when stratified by prediction set size (coverage
is only shown if there was at least one prediction set at the desired $\alpha$
for a particular method).

We also compare coverage and distribution stratified by set size in Figure 6.
Stratifying by set size reveals that most predictions with Ordinal APS and
Naive LAC contain only one or two grading predictions while Ordinal CDF mainly
predicts a set of all four possible gradings (which always trivially satisfies
coverage).

Lastly, we compare coverage and set size at the disc level in Table 1 at
$\alpha=0.1$. We find that coverage was inversely correlated to the prevalence
of severe stenosis, which is most often found in the lower lumbar disc levels.

Table 1: Coverage and set size stratified by intervertebral disc level at $\alpha=0.1$. disc level | average severity | method | coverage | set size
---|---|---|---|---
T12-L1 | 0.04 | Ordinal CDF | $99.7\%\pm 0.3\%$ | $3.97\pm 0.02$
Naive LAC | $98.6\%\pm 0.4\%$ | $1.04\pm 0.02$
Ordinal APS | $98.6\%\pm 0.3\%$ | $1.04\pm 0.01$
L1-L2 | 0.18 | Ordinal CDF | $97.5\%\pm 1.1\%$ | $3.84\pm 0.06$
Naive LAC | $95.4\%\pm 1.0\%$ | $1.12\pm 0.04$
Ordinal APS | $95.4\%\pm 1.0\%$ | $1.12\pm 0.04$
L2-L3 | 0.48 | Ordinal CDF | $90.4\%\pm 2.1\%$ | $3.50\pm 0.11$
Naive LAC | $90.0\%\pm 3.0\%$ | $1.25\pm 0.09$
Ordinal APS | $90.1\%\pm 3.0\%$ | $1.26\pm 0.09$
L3-L4 | 0.75 | Ordinal CDF | $84.3\%\pm 4.0\%$ | $3.17\pm 0.19$
Naive LAC | $86.0\%\pm 4.0\%$ | $1.42\pm 0.12$
Ordinal APS | $85.7\%\pm 4.2\%$ | $1.42\pm 0.13$
L4-L5 | 1.06 | Ordinal CDF | $81.7\%\pm 4.2\%$ | $2.86\pm 0.21$
Naive LAC | $83.2\%\pm 4.6\%$ | $1.48\pm 0.15$
Ordinal APS | $83.2\%\pm 4.8\%$ | $1.48\pm 0.15$
L5-S1 | 0.71 | Ordinal CDF | $84.4\%\pm 4.1$ | $3.20\pm 0.17$
Naive LAC | $85.3\%\pm 3.2$ | $1.39\pm 0.12$
Ordinal APS | $85.7\%\pm 3.0$ | $1.41\pm 0.11$
Total | 0.54 | Ordinal CDF | $90.0\%\pm 2.5\%$ | $3.44\pm 0.12$
Naive LAC | $90.0\%\pm 2.6\%$ | $1.28\pm 0.09$
Ordinal APS | $90.1\%\pm 2.7\%$ | $1.28\pm 0.09$

Overall, we conclude that Ordinal APS performs similarly to LAC in both
coverage and set size, and both Ordinal APS and Naive LAC outperform Ordinal
CDF. The similarities between LAC and Ordinal APS are notable — they almost
always result in the same sets, although the algorithms are quite different.
This is unsurprising given that in our setting $|\mathcal{Y}|=4$ and the
model’s accuracy is high, so bimodal softmax scores almost never happen.
Therefore LAC and Ordinal APS do the same thing. This observation does not
generalize; other ordinal prediction problems with more categories will have
bimodal distributions and thus LAC and Ordinal APS will differ.

###  3.2 Clinical Review of High Uncertainty Predictions

To investigate the clinical utility of Ordinal APS to enhance AI-augmented
workflows, we evaluate one possible clinical integration use case: flagging
low confidence predictions (i.e. ones with a large set size). The
radiologist’s performance and user experience of the model may be improved by
raising their awareness of those instances in which the model performance may
be degraded by scan anomalies or when uncertainty quantification is very high,
excluding such potentially poor quality results from their review
responsibilities altogether. We define an uncertainty score for each patient
by taking the average set size for all disc levels and grading tasks. An
neuroradiologist with $>20$ years of experience determined what constituted a
“significant imaging anomaly” within the context of spinal stenosis
interpretation.

Figure 7: Three anomalies found in high uncertainty predictions; the anomalous
areas are boxed in red.

As a statistical validation of these results, we examined the report of 70
cases with the highest uncertainty and found 17 such anomalies: 11 cases with
artifacts from metallic orthopedic hardware, four cases with motion artifacts,
one case with a large tumor occupying the spinal canal, and one case with a
severe congenital abnormality (achondroplastic dwarfism). In contrast, a
random sample of 70 cases from the dataset only demonstrated five cases with
significant anomalies which were all orthopedic hardware artifacts. This
difference is significant with $p<0.05$ by Fisher’s exact test, and
qualitatively the abnormalities found in the filtered samples were more
extreme.

Our manual review of high uncertainty cases shows promise for improving the
clinician experience with AI-assisted medical software tools using
distribution-free uncertainty quantification. Rather than presenting all AI
predictions as equally trustworthy, cases with higher uncertainty can be
flagged to prompt additional scrutiny or hidden from the user altogether to
maximize efficiency. While prospective evaluation of this use of conformal
prediction in more clinically realistic settings will be needed to validate
general feasibility and utility, our preliminary experiments are a step
towards demonstrating the clinical applicability of conformal prediction.

##  4 Related Work

#### Conformal Prediction and Distribution-Free Uncertainty Quantification

Conformal prediction is a flexible technique for generating prediction
intervals from arbitrary models. It was first developed by Vladimir Vovk and
collaborators in the late 1990s [32, 31, 17, 18, 16, 26]. We build most
directly on the work of Yaniv Romano and collaborators, who developed the
Adaptive Prediction Sets method studied in [25, 6] and the histogram binning
method in [27]. The latter work is the most relevant, and it proposes an
algorithm very similar to Algorithm 1 in a continuous setting with histogram
regression. Our work also relies on the nested set outlook on conformal
prediction [13]. We also build directly on existing work involving
distribution-free risk-controlling prediction sets and Learn then Test [8, 3].
The LAC baseline is taken from [26], and the ordinal CDF baseline is similar
to the softmax method in [5], which is in turn motivated by [15, 24]. A gentle
introduction to these topics and their history is available in [2], or
alternatively, in [28].

#### Uncertainty Quantification for Critical Applications

Recently, uncertainty quantification has been promoted to facilitate
trustworthiness and transparency in black-box algorithms, such as deep
learning, for critical decision-making [10]. In particular, conformal
prediction methods have been applied to a wide range of safety-critical
applications – from reducing false alarms in the detection of sepsis risk [29]
to end-to-end autonomous driving control [22]. Distribution-free uncertainty
quantification techniques such as conformal prediction sets have emerged as an
essential tool for rigorous statistical guarantees in medical decision-making
[9, 20, 4, 30, 11, 19].

##  5 Conclusion

We show how conformal prediction sets can complement existing AI systems
without further modification to the model to provide distribution-free
reliability guarantees. We demonstrate its clinically utility in the
application of flagging high uncertainty cases in automated stenosis severity
grading for followup review. We hope this work promotes further studies on the
trustworthiness and usability of uncertainty-aware machine learning systems
for clinical applications.

## Acknowledgements

SP receives research support from GE Healthcare. AA is funded by the NSFGRFP
and a Berkeley Fellowship. We thank Stephen Bates for many helpful
conversations.

## References

  * [1] Allen, B., Agarwal, S., Coombs, L., Wald, C., Dreyer, K.: 2020 acr data science institute artificial intelligence survey. Journal of the American College of Radiology 18(8), 1153–1159 (2021)
  * [2] Angelopoulos, A.N., Bates, S.: A gentle introduction to conformal prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.07511 (2021)
  * [3] Angelopoulos, A.N., Bates, S., Candès, E.J., Jordan, M.I., Lei, L.: Learn then test: Calibrating predictive algorithms to achieve risk control. arXiv preprint arXiv:2110.01052 (2021)
  * [4] Angelopoulos, A.N., Bates, S., Zrnic, T., Jordan, M.I.: Private prediction sets. arXiv preprint arXiv:2102.06202 (2021)
  * [5] Angelopoulos, A.N., Kohli, A.P., Bates, S., Jordan, M.I., Malik, J., Alshaabi, T., Upadhyayula, S., Romano, Y.: Image-to-image regression with distribution-free uncertainty quantification and applications in imaging. arXiv preprint arXiv:2202.05265 (2022)
  * [6] Angelopoulos, A.N., Bates, S., Malik, J., Jordan, M.I.: Uncertainty sets for image classifiers using conformal prediction. In: International Conference on Learning Representations (ICLR) (2021)
  * [7] Barber, R., Candès, E., Ramdas, A., Tibshirani, R.: The limits of distribution-free conditional predictive inference. Information and Inference 10(2), 455–482 (08 2021)
  * [8] Bates, S., Angelopoulos, A., Lei, L., Malik, J., Jordan, M.: Distribution-free, risk-controlling prediction sets. Journal of the Association for Computing Machinery 68(6) (9 2021)
  * [9] Begoli, E., Bhattacharya, T., Kusnezov, D.: The need for uncertainty quantification in machine-assisted medical decision making. Nature (01 2019)
  * [10] Bhatt, U., Antorán, J., Zhang, Y., Liao, Q.V., Sattigeri, P., Fogliato, R., Melançon, G., Krishnan, R., Stanley, J., Tickoo, O., Nachman, L., Chunara, R., Srikumar, M., Weller, A., Xiang, A.: Uncertainty as a Form of Transparency: Measuring, Communicating, and Using Uncertainty, p. 401–413. Association for Computing Machinery, New York, NY, USA (2021)
  * [11] Fannjiang, C., Bates, S., Angelopoulos, A., Listgarten, J., Jordan, M.I.: Conformal prediction for the design problem. arXiv preprint arXiv:2202.03613 (2022)
  * [12] Fu, M., Buerba, R., III, W., Blizzard, D., Lischuk, A., Haims, A., Grauer, J.: Inter-rater and intra-rater agreement of magnetic resonance imaging findings in the lumbar spine: Significant variability across degenerative conditions. The Spine Journal 14 (10 2014). https://doi.org/10.1016/j.spinee.2014.03.010
  * [13] Gupta, C., Kuchibhotla, A.K., Ramdas, A.: Nested conformal prediction and quantile out-of-bag ensemble methods. Pattern Recognition p. 108496 (2021)
  * [14] Hendrycks, D., Zhao, K., Basart, S., Steinhardt, J., Song, D.: Natural adversarial examples. CVPR (2021)
  * [15] Koenker, R., Bassett Jr, G.: Regression quantiles. Econometrica: Journal of the Econometric Society 46(1), 33–50 (1978)
  * [16] Lei, J.: Classification with confidence. Biometrika 101(4), 755–769 (10 2014)
  * [17] Lei, J., Rinaldo, A., Wasserman, L.: A conformal prediction approach to explore functional data. Annals of Mathematics and Artificial Intelligence 74, 29–43 (2015)
  * [18] Lei, J., Robins, J., Wasserman, L.: Distribution-free prediction sets. Journal of the American Statistical Association 108(501), 278–287 (2013)
  * [19] Lu, C., Chang, K., Singh, P., Kalpathy-Cramer, J.: Three applications of conformal prediction for rating breast density in mammography (2022). https://doi.org/10.48550/ARXIV.2206.12008, https://arxiv.org/abs/2206.12008
  * [20] Lu, C., Lemay, A., Chang, K., Hoebel, K., Kalpathy-Cramer, J.: Fair conformal predictors for applications in medical imaging (2022)
  * [21] Lu, J., Pedemonte, S., Bizzo, B., Doyle, S., Andriole, K.P., Michalski, M.H., Gonzalez, R.G., Pomerantz, S.R.: Deepspine: Automated lumbar vertebral segmentation, disc-level designation, and spinal stenosis grading using deep learning. CoRR abs/1807.10215 (2018)
  * [22] Michelmore, R., Wicker, M., Laurenti, L., Cardelli, L., Gal, Y., Kwiatkowska, M.: Uncertainty quantification with statistical guarantees in end-to-end autonomous driving control. CoRR abs/1909.09884 (2019)
  * [23] Minderer, M., Djolonga, J., Romijnders, R., Hubis, F.A., Zhai, X., Houlsby, N., Tran, D., Lucic, M.: Revisiting the calibration of modern neural networks. arXiv preprint arXiv:2106.07998 (2021)
  * [24] Romano, Y., Patterson, E., Candès, E.: Conformalized quantile regression. In: Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, vol. 32, pp. 3543–3553. NIPS (2019)
  * [25] Romano, Y., Sesia, M., Candès, E.: Classification with valid and adaptive coverage. In: Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M.F., Lin, H. (eds.) Advances in Neural Information Processing Systems. vol. 33, pp. 3581–3591. Curran Associates, Inc. (2020)
  * [26] Sadinle, M., Lei, J., Wasserman, L.: Least ambiguous set-valued classifiers with bounded error levels. Journal of the American Statistical Association 114, 223 – 234 (2019)
  * [27] Sesia, M., Romano, Y.: Conformal prediction using conditional histograms. Advances in Neural Information Processing Systems 34 (2021)
  * [28] Shafer, G., Vovk, V.: A tutorial on conformal prediction. Journal of Machine Learning Research 9(Mar), 371–421 (2008)
  * [29] Shashikumar, S., Wardi, G., Malhotra, A., Nemati, S.: Artificial intelligence sepsis prediction algorithm learns to say “I don’t know”. npj Digital Medicine 4 (12 2021)
  * [30] Vazquez, J., Facelli, J.: Conformal prediction in clinical medical sciences. Journal of Healthcare Informatics Research (01 2022)
  * [31] Vovk, V., Gammerman, A., Shafer, G.: Algorithmic Learning in a Random World. Springer, New York, NY, USA (2005)
  * [32] Vovk, V., Gammerman, A., Saunders, C.: Machine-learning applications of algorithmic randomness. In: International Conference on Machine Learning. pp. 444–453 (1999)

##  Appendix A Formal Theorem Statements and Proofs

######  Theorem 1 (Formal conformal coverage guarantee).

Let $(X_{1},Y_{1})$, $(X_{2},Y_{2})$, …, $(X_{n},Y_{n})$ and $(X_{\rm
test},Y_{\rm test})$ be drawn independently and identically distributed from
distribution $\mathbb{P}$, and let $\mathcal{C}_{\lambda}$ be any sequence of
sets nested in $\lambda$ such that
$\underset{\lambda\to\infty}{\lim}\mathcal{C}_{\lambda}=\mathcal{Y}$. Finally,
let $\hat{\lambda}=\mathcal{A}(\mathcal{C}_{\lambda},\;\alpha)$. Then
$\mathcal{C}_{\hat{\lambda}}$ satisfies (3).

###### Theorem 1.

This is a standard result in the conformal literature [32]. A stripped down
proof appears in [2]. The nested set version appears in [13]. ∎

The informal version of Theorem 1 in the main text is simply a corollary of
the above when applied to the specific sequence of nested sets
$\mathcal{T}_{\lambda}$.

###### Corollary 1.

Pick $\lambda^{(1)}\leq\lambda^{(2)}$. It is clear that
$\hat{f}(j|x)\geq\lambda^{(2)}$ implies $\hat{f}(j|x)\geq\lambda^{(1)}$.
Therefore
$\mathcal{T}_{\lambda^{(2)}}(x)\subseteq\mathcal{T}_{\lambda^{(1)}}(x)$.
Applying Theorem 1 completes the proof. ∎

######  Corollary A.1.

In the setting of Theorem 1, let $\tilde{\mathcal{T}}_{\lambda}$ be the
sequence of nested sets computed using Algorithm 1, and let
$\hat{\lambda}=\mathcal{A}\Big{(}\tilde{\mathcal{T}}_{\lambda},\;\alpha\Big{)}$.
Then $\tilde{\mathcal{T}}_{\hat{\lambda}}$ satisfies 3.

###### Corollary A.1.

Pick $\lambda^{(1)}\leq\lambda^{(2)}$. Examining Algorithm 1,
$q\leq\lambda^{(1)}$ implies $q\leq\lambda^{(2)}$ . Therefore
$\tilde{\mathcal{T}}_{\lambda^{(1)}}(x)\subseteq\tilde{\mathcal{T}}_{\lambda^{(2)}}(x)$.
Applying Theorem 1 completes the proof. ∎