\mathbf{V}_{\perp}^{*}\end{bmatrix},$

e.g., see [Van13, Proposition 2.1]. We thus observe that the weight operator
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ from Eq. 19 satisfies

$\begin{split}&W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{M})\\\
&=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\left(\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})\circ\left(\begin{bmatrix}\mathbf{U}^{*}\\\
\mathbf{U}_{\perp}^{*}\end{bmatrix}\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\begin{bmatrix}\mathbf{M}_{1}&\mathbf{M}_{2}\\\
\mathbf{M}_{3}&0\end{bmatrix}\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}\begin{bmatrix}\mathbf{V}&\mathbf{V}_{\perp}\end{bmatrix}\right)\right)\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}\\\
&=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\left(\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})\circ\begin{bmatrix}\mathbf{M}_{1}&\mathbf{M}_{2}\\\
\mathbf{M}_{3}&0\end{bmatrix}\right)\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}\\\
&=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\begin{bmatrix}\mathbf{H}_{1}^{(k)}\circ\mathbf{M}_{1}&\mathbf{H}_{2}^{(k)}\circ\mathbf{M}_{2}\\\
\mathbf{H}_{3}^{(k)}\circ\mathbf{M}_{3}&0\end{bmatrix}\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}\in T_{\mathbf{U},\mathbf{V}}.\end{split}$

Similarly, if $\mathbf{M}\in T_{\mathbf{U},\mathbf{V}}^{\perp}$, there exists
$\mathbf{M}_{4}\in\mathbb{R}^{(n_{1}-r_{k})\times(n_{2}-r_{k})}$ such that
$\mathbf{M}=\mathbf{U}_{\perp}\mathbf{M}_{4}(\mathbf{V}_{\perp})^{*}$ and

$\begin{split}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{M})&=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\left(\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})\circ\begin{bmatrix}0&0\\\
0&\mathbf{M}_{4}\end{bmatrix}\right)\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}\\\
&=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\begin{bmatrix}0&0\\\
0&\mathbf{M}_{4}/\varepsilon_{k}^{2}\end{bmatrix}\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}\in
T_{\mathbf{U},\mathbf{V}}^{\perp}.\end{split}$

######  Proof D.4 (of Lemma C.7).

Let $\boldsymbol{\Xi}\in\mathbb{R}^{n_{1}\times n_{2}}$ be arbitrary. We start
with some simple but technical observations. First note that by Lemma C.1, if
$T_{k}=T_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}},S}$,

(45)
$\displaystyle\begin{split}\mathbb{P}_{T_{k}^{\perp}}\boldsymbol{\Xi}&=(\mathbf{Id}-\mathbb{P}_{T_{k}})(\boldsymbol{\Xi})\\\
&=\boldsymbol{\Xi}-\Big{(}\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*}\boldsymbol{\Xi}+\mathbb{P}_{S}\boldsymbol{\Xi}\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*}\boldsymbol{\Xi}\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}\Big{)}\\\
&=(\mathbf{Id}-\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}})\boldsymbol{\Xi}+(\mathbf{Id}-\mathbb{P}_{S})\boldsymbol{\Xi}\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}\\\
&=\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}\boldsymbol{\Xi}+\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*},\end{split}$

with

(46)
$\displaystyle\Big{\langle}\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}\boldsymbol{\Xi},\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}\Big{\rangle}=\Big{\langle}\widetilde{\mathbf{U}}_{\perp}\widetilde{\mathbf{U}}_{\perp}^{*}\boldsymbol{\Xi}\widetilde{\mathbf{V}}_{\perp}\widetilde{\mathbf{V}}_{\perp}^{*},\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}\Big{\rangle}=0,$

where we used that
$\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}\boldsymbol{\Xi}=(\mathbf{Id}-\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}})\boldsymbol{\Xi}=\widetilde{\mathbf{U}}_{\perp}\widetilde{\mathbf{U}}_{\perp}^{*}\boldsymbol{\Xi}\widetilde{\mathbf{V}}_{\perp}\widetilde{\mathbf{V}}_{\perp}^{*}$,
for $\widetilde{\mathbf{U}}_{\perp}\in\mathbb{R}^{n_{1}\times(n_{1}-r)}$ and
$\widetilde{\mathbf{V}}_{\perp}\in\mathbb{R}^{n_{2}\times(n_{2}-r)}$ being the
complementary orthonormal bases of $\widetilde{\mathbf{U}}$ and
$\widetilde{\mathbf{V}}$, and that
$\widetilde{\mathbf{V}}_{\perp}^{*}\widetilde{\mathbf{V}}=\boldsymbol{0}$.

Second, let now $\mathbf{U}\in\mathbb{R}^{n_{1}\times r}$ and
$\mathbf{V}\in\mathbb{R}^{n_{2}\times r}$ be matrices with $r$ leading left
and right singular vectors of $\mathbf{X}^{(k)}$ in their columns which
coincide with the matrices $\mathbf{U}$ and $\mathbf{V}$ from Definition 2.1
in their first $r$ columns. Then it follows from Eqs. 45 and 46 that

(47)
$\displaystyle\begin{split}&\left\|\mathbb{P}_{T_{k}^{\perp}}(\boldsymbol{\Xi})\right\|_{F}^{2}\\\
&=\left\|\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}\boldsymbol{\Xi}\right\|_{F}^{2}+\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}\right\|_{F}^{2}\\\
&=\left\|\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}+(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}-\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp})\boldsymbol{\Xi}\right\|_{F}^{2}+\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}(\mathbf{V}\mathbf{V}^{*}+(\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}\mathbf{V}^{*}))\right\|_{F}^{2}\\\
&\leq
2\left(\left\|\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}\right\|_{F}^{2}+\left\|(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}-\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp})\boldsymbol{\Xi}\right\|_{F}^{2}+\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\mathbf{V}\mathbf{V}^{*}\right\|_{F}^{2}+\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}(\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}\mathbf{V}^{*})\right\|_{F}^{2}\right).\end{split}$

By an argument analogous to Eq. 46, we observe that

$\left\|\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}\right\|_{F}^{2}+\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\mathbf{V}\mathbf{V}^{*}\right\|_{F}^{2}=\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}+\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\mathbf{V}\mathbf{V}^{*},\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}+\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\mathbf{V}\mathbf{V}^{*}\rangle=\langle\widetilde{\boldsymbol{\Xi}},\widetilde{\boldsymbol{\Xi}}\rangle,$

where
$\widetilde{\boldsymbol{\Xi}}=\mathbb{P}_{\mathcal{M}}(\boldsymbol{\Xi})$ is
an element of the subspace
$\mathcal{M}=\mathcal{M}_{1}\oplus\mathcal{M}_{2}\subset\mathbb{R}^{n_{1}\times
n_{2}}$ that is the direct sum of the subspaces
$\mathcal{M}_{1}:=\\{\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\mathbf{Z}:\mathbf{Z}\in\mathbb{R}^{n_{1}\times
n_{2}}\\}$ and
$\mathcal{M}_{2}:=\\{\mathbb{P}_{S^{c}}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}:\mathbf{Z}\in\mathbb{R}^{n_{1}\times
n_{2}}\\}$.

Let now $W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ be the rank
promoting part of the weight operator from Eq. 9 and
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$ be the row-sparsity
promoting part from Eq. 10. Note that the restriction of

$\bar{W}:=\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}+\mathbb{P}_{S^{c}}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}$

to $\mathcal{M}$ is invertible as its first summand is invertible on
$\mathcal{M}_{1}=T_{\mathbf{U},\mathbf{V}}^{\perp}$, its second summand is
invertible on $\mathcal{M}_{2}$ (recall that the weight operators are positive
definite), and $\mathcal{M}_{1}\perp\mathcal{M}_{2}$. Therefore it holds that

$\begin{split}&\left\|\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}\right\|_{F}^{2}+\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\mathbf{V}\mathbf{V}^{*}\right\|_{F}^{2}\\\
&=\langle\widetilde{\boldsymbol{\Xi}},\widetilde{\boldsymbol{\Xi}}\rangle=\Big{\langle}\bar{W}_{|\mathcal{M}}^{1/2}\widetilde{\boldsymbol{\Xi}},\bar{W}_{|\mathcal{M}}^{-1}\bar{W}_{|\mathcal{M}}^{1/2}\widetilde{\boldsymbol{\Xi}}\Big{\rangle}\\\
&\leq\sigma_{1}\left(\bar{W}_{|\mathcal{M}}^{-1}\right)\Big{\langle}\widetilde{\boldsymbol{\Xi}},\bar{W}_{|\mathcal{M}}\widetilde{\boldsymbol{\Xi}}\Big{\rangle}=\frac{1}{\sigma_{\operatorname{min}}\left(\bar{W}_{|\mathcal{M}}\right)}\Big{\langle}\widetilde{\boldsymbol{\Xi}},\bar{W}\widetilde{\boldsymbol{\Xi}}\Big{\rangle}\\\
&\leq\frac{1}{\sigma_{\operatorname{min}}\left(\left(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\right)_{|\mathcal{M}}\right)+\sigma_{\operatorname{min}}\left(\left(\mathbb{P}_{S^{c}}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\right)_{|\mathcal{M}}\right)}\Big{\langle}\widetilde{\boldsymbol{\Xi}},\bar{W}\widetilde{\boldsymbol{\Xi}}\Big{\rangle}\\\
&\leq\frac{1}{\frac{\varepsilon_{k}^{2}}{\sigma_{r+1}^{2}(\mathbf{X}^{(k)})}+\frac{\delta_{k}^{2}}{\rho_{s+1}^{2}(\mathbf{X}^{(k)})}}\Big{\langle}\boldsymbol{\Xi},\bar{W}\boldsymbol{\Xi}\Big{\rangle}.\end{split}$

In the first inequality, we used that $\bar{W}_{|\mathcal{M}}$ is positive
definite. In the second inequality, we used that
$\sigma_{\min}(A+B)\geq\sigma_{\min}(A)+\sigma_{\min}(B)$, for any positive
semidefinite operators $A$ and $B$, and and in the third inequality that
$\langle\widetilde{\boldsymbol{\Xi}},\bar{W}\widetilde{\boldsymbol{\Xi}}\rangle\leq\langle\boldsymbol{\Xi},\bar{W}\boldsymbol{\Xi}\rangle$.
The latter observation can be deduced as follows: Note that, by the self-
adjointness of $\bar{W}$,

$\displaystyle\langle\boldsymbol{\Xi},\bar{W}\boldsymbol{\Xi}\rangle$
$\displaystyle=\langle\mathbb{P}_{\mathcal{M}}(\boldsymbol{\Xi}),\bar{W}\mathbb{P}_{\mathcal{M}}(\boldsymbol{\Xi})\rangle+\langle\mathbb{P}_{\mathcal{M}}(\boldsymbol{\Xi}),\bar{W}\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi})\rangle+\langle\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi}),\bar{W}\mathbb{P}_{\mathcal{M}}(\boldsymbol{\Xi})\rangle+\langle\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi}),\bar{W}\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi})\rangle$
$\displaystyle=\langle\widetilde{\boldsymbol{\Xi}},\bar{W}\widetilde{\boldsymbol{\Xi}}\rangle+2\langle\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi}),\bar{W}\mathbb{P}_{\mathcal{M}}(\boldsymbol{\Xi})\rangle+\langle\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi}),\bar{W}\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi})\rangle.$

Since $\bar{W}$ is positive semi-definite (due to the fact that both
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ and
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$ are positive
definite), all that remains is to argue that the mixed term on the right-hand
side vanishes. To this end, note that
$\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\mathbf{Z}=\mathbb{P}_{S^{c}}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}=0$,
for any $\mathbf{Z}\in\mathcal{M}_{\perp}$ and compute

$\displaystyle\langle\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi}),\bar{W}\mathbb{P}_{\mathcal{M}}(\boldsymbol{\Xi})\rangle$
$\displaystyle=\langle\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi}),\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\mathbb{P}_{\mathcal{M}_{1}}(\boldsymbol{\Xi}))))\rangle+\langle\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi}),\mathbb{P}_{S^{c}}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\mathbb{P}_{\mathcal{M}_{2}}(\boldsymbol{\Xi})\rangle$
$\displaystyle=\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi})),W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\mathbb{P}_{\mathcal{M}_{1}}(\boldsymbol{\Xi})))\rangle+\langle\mathbb{P}_{S^{c}}\mathbb{P}_{\mathcal{M}}^{\perp}(\boldsymbol{\Xi})\mathbf{V}\mathbf{V}^{*},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\rangle$
$\displaystyle=0.$

We can now continue by estimating

(48)
$\displaystyle\begin{split}\Big{\langle}\boldsymbol{\Xi},\bar{W}\boldsymbol{\Xi}\Big{\rangle}&=\Big{\langle}\boldsymbol{\Xi},\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}\Big{\rangle}+\Big{\langle}\boldsymbol{\Xi},\mathbb{P}_{S^{c}}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\Big{\rangle}\\\
&=\Big{\langle}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\boldsymbol{\Xi}),W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\boldsymbol{\Xi})\Big{\rangle}+\Big{\langle}\mathbb{P}_{S^{c}}\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\Big{\rangle}\\\
&\leq\Big{\langle}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\boldsymbol{\Xi}\Big{\rangle}+\Big{\langle}\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\boldsymbol{\Xi}\Big{\rangle}\\\
&=\Big{\langle}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}\boldsymbol{\Xi}\Big{\rangle},\end{split}$

using the positive semidefiniteness of
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ and
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$ in the last
inequality. To be precise, the last inequality can be argued as follows: Due
to complimentary supports $S$ and $S^{c}$, we see that

(49)
$\displaystyle\begin{split}\langle\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\boldsymbol{\Xi}\rangle&=\langle\mathbb{P}_{S}\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S}\boldsymbol{\Xi}\rangle+\langle\mathbb{P}_{S^{c}}\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\rangle\\\
&+\underbrace{\langle\mathbb{P}_{S}\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\rangle}_{=0}+\underbrace{\langle\mathbb{P}_{S^{c}}\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S}\boldsymbol{\Xi}\rangle}_{=0}\\\
&\geq\langle\mathbb{P}_{S^{c}}\boldsymbol{\Xi},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}\rangle.\end{split}$

Similarly, we note that $W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$
acts diagonally on $T_{\mathbf{U},\mathbf{V}}$ and
$T_{\mathbf{U},\mathbf{V}}^{\perp}$. Indeed, we have by Lemma D.2 that if
$\mathbf{M}\in T_{\mathbf{U},\mathbf{V}}$, then
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{M})\in
T_{\mathbf{U},\mathbf{V}}$ and if $\mathbf{M}\in
T_{\mathbf{U},\mathbf{V}}^{\perp}$, then
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{M})\in
T_{\mathbf{U},\mathbf{V}}^{\perp}$, which implies

$\displaystyle\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi})\rangle=0\quad\text{and}\quad\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}\boldsymbol{\Xi})\rangle=0$

due to the orthogonality of elements in $T_{\mathbf{U},\mathbf{V}}^{\perp}$
and $T_{\mathbf{U},\mathbf{V}}$, respectively, and therefore it follows from
$\boldsymbol{\Xi}=T_{\mathbf{U},\mathbf{V}}(\boldsymbol{\Xi})+T_{\mathbf{U},\mathbf{V}}^{\perp}(\boldsymbol{\Xi})$
that

(50)
$\displaystyle\begin{split}\langle\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\boldsymbol{\Xi})\rangle&=\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}\boldsymbol{\Xi})\rangle+\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi})\rangle\\\
&+\underbrace{\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi})\rangle}_{=0}+\underbrace{\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}\boldsymbol{\Xi})\rangle}_{=0}\\\
&\geq\langle\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi})\rangle.\end{split}$

Combining the previous estimates with Eq. 47 and noticing that

$\frac{1}{\frac{\varepsilon_{k}^{2}}{\sigma_{r+1}^{2}(\mathbf{X}^{(k)})}+\frac{\delta_{k}^{2}}{\rho_{s+1}^{2}(\mathbf{X}^{(k)})}}\leq\min\left\\{\frac{\sigma_{r+1}^{2}(\mathbf{X}^{(k)})}{\varepsilon_{k}^{2}},\frac{\rho_{s+1}^{2}(\mathbf{X}^{(k)})}{\delta_{k}^{2}}\right\\},$

we obtain

(51)
$\begin{split}\left\|\mathbb{P}_{T_{k}^{\perp}}(\boldsymbol{\Xi})\right\|_{F}^{2}&\leq
2\left(\min\left\\{\frac{\sigma_{r+1}^{2}(\mathbf{X}^{(k)})}{\varepsilon_{k}^{2}},\frac{\rho_{s+1}^{2}(\mathbf{X}^{(k)})}{\delta_{k}^{2}}\right\\}\Big{\langle}\boldsymbol{\Xi},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}\boldsymbol{\Xi}\Big{\rangle}\right)\\\
&+2\left\|(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}-\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp})\boldsymbol{\Xi}\right\|_{F}^{2}+2\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}(\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}\mathbf{V}^{*})\right\|_{F}^{2}.\end{split}$

Next, we control the last two summands in Eq. 51 using matrix perturbation
results. Recall that $\mathbf{U}_{\star}\in\mathbb{R}^{n_{1}\times r}$ and
$\mathbf{V}_{\star}\in\mathbb{R}^{n_{2}\times r}$ are the singular vector
matrices of the reduced singular value decomposition of $\mathbf{X}_{\star}$.
First observe that

$\begin{split}&\left\|(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}-\mathbb{P}_{\mathbf{U}_{\star},\mathbf{V}_{\star}})\boldsymbol{\Xi}\right\|_{F}\\\
&=\left\|(\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*}-\mathbf{U}_{\star}\mathbf{U}_{\star}^{*})\boldsymbol{\Xi}(\mathbf{Id}-\mathbf{V}_{\star}\mathbf{V}_{\star}^{*})+(\mathbf{Id}-\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*})\boldsymbol{\Xi}(\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}_{\star}\mathbf{V}_{\star}^{*})\right\|_{F}\\\
&\leq\left\|\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*}-\mathbf{U}_{\star}\mathbf{U}_{\star}^{*}\right\|\left\|\boldsymbol{\Xi}\right\|_{F}\left\|\mathbf{Id}-\mathbf{V}_{\star}\mathbf{V}_{\star}^{*}\right\|+\left\|\mathbf{Id}-\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*}\right\|\left\|\boldsymbol{\Xi}\right\|_{F}\left\|\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}_{\star}\mathbf{V}_{\star}^{*}\right\|\\\
&\leq\left(\left\|\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*}-\mathbf{U}_{\star}\mathbf{U}_{\star}^{*}\right\|+\left\|\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}_{\star}\mathbf{V}_{\star}^{*}\right\|\right)\left\|\boldsymbol{\Xi}\right\|_{F}.\end{split}$

Now note that, by [CZ18, Lemma 1] and [LW20, Theorem 3.5], we obtain

(52)
$\begin{split}\left\|\widetilde{\mathbf{U}}\widetilde{\mathbf{U}}^{*}-\mathbf{U}_{\star}\mathbf{U}_{\star}^{*}\right\|+\left\|\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}_{\star}\mathbf{V}_{\star}^{*}\right\|&\leq
2\left(\|\mathbf{U}_{\star}\widetilde{\mathbf{U}}_{\perp}^{*}\|+\|\mathbf{V}_{\star}\widetilde{\mathbf{V}}_{\perp}^{*}\|\right)\leq
4\frac{\|\mathsf{H}_{s}(\mathbf{X}^{(k)})-\mathbf{X}_{\star}\|}{\sigma_{r}(\mathbf{X}_{\star})-\sigma_{r+1}(\mathsf{H}_{s}(\mathbf{X}^{(k)}))}\\\
&\leq
4\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{\sigma_{r}(\mathbf{X}_{\star})-\sigma_{r+1}(\mathbf{X}^{(k)})}\leq
4\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{(1-1/48)\sigma_{r}(\mathbf{X}_{\star})},\end{split}$

where we used some small observations in the last two inequalities: First,
$\sigma_{r+1}(\mathsf{H}_{s}(\mathbf{X}^{(k)}))\leq\sigma_{r+1}(\mathbf{X}^{(k)})$,
which follows from the rectangular Cauchy interlacing theorem [Dax10, Theorem
23]. Second, according to Lemma C.2 and Eq. 26, the row-support $S$ of
$\mathsf{H}_{s}(\mathbf{X}^{(k)})$ coincides with the row-support
$S_{\star}=\\{i\in[n_{1}]:\|(\mathbf{X}_{\star})_{i,:}\|_{2}\neq 0\\}$ of
$\mathbf{X}_{\star}$ and hence
$\mathsf{H}_{s}(\mathbf{X}^{(k)})-\mathbf{X}_{\star}$ is a submatrix of
$\mathbf{X}^{(k)}-\mathbf{X}_{\star}$. Finally,
$\sigma_{r+1}(\mathbf{X}^{(k)})=\|\mathsf{T}_{r}(\mathbf{X}^{(k)})-\mathbf{X}^{(k)}\|\leq\|\mathbf{X}_{\star}-\mathbf{X}^{(k)}\|\leq\frac{1}{48}\sigma_{r}(\mathbf{X}_{\star})$
due to Eq. 26. Consequently,

$\displaystyle\left\|(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}-\mathbb{P}_{\mathbf{U}_{\star},\mathbf{V}_{\star}})\boldsymbol{\Xi}\right\|_{F}\leq
4\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{(1-1/48)\sigma_{r}(\mathbf{X}_{\star})}\left\|\boldsymbol{\Xi}\right\|_{F}$

and, by a similar argument,

$\displaystyle\left\|(\mathbb{P}_{\mathbf{U}_{\star},\mathbf{V}_{\star}}-\mathbb{P}_{\mathbf{U},\mathbf{V}})\boldsymbol{\Xi}\right\|_{F}\leq
4\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{(1-1/48)\sigma_{r}(\mathbf{X}_{\star})}\left\|\boldsymbol{\Xi}\right\|_{F},$

such that it follows from
$\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}=\mathbf{Id}-\mathbb{P}_{\mathbf{U},\mathbf{V}}$
and
$\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}=\mathbf{Id}-\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}$
that

(53)
$\begin{split}\left\|(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}^{\perp}-\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp})\boldsymbol{\Xi}\right\|_{F}&=\left\|(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}-\mathbb{P}_{\mathbf{U},\mathbf{V}})\boldsymbol{\Xi}\right\|_{F}\leq\left\|(\mathbb{P}_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}}}-\mathbb{P}_{\mathbf{U}_{\star},\mathbf{V}_{\star}})\boldsymbol{\Xi}\right\|_{F}+\Big{\|}(\mathbb{P}_{\mathbf{U}_{\star},\mathbf{V}_{\star}}-\mathbb{P}_{\mathbf{U},\mathbf{V}})\boldsymbol{\Xi}\Big{\|}_{F}\\\
&\leq
8\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{(1-1/48)\sigma_{r}(\mathbf{X}_{\star})}\left\|\boldsymbol{\Xi}\right\|_{F}.\end{split}$

To estimate the fourth summand in Eq. 51, we argue analogously that

(54)
$\begin{split}\left\|\mathbb{P}_{S^{c}}\boldsymbol{\Xi}(\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}\mathbf{V}^{*})\right\|_{F}&\leq\left\|\boldsymbol{\Xi}(\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}\mathbf{V}^{*})\right\|_{F}\leq\left(\|\widetilde{\mathbf{V}}\widetilde{\mathbf{V}}^{*}-\mathbf{V}_{\star}\mathbf{V}_{\star}^{*}\|+\|\mathbf{V}_{\star}\mathbf{V}_{\star}^{*}-\mathbf{V}\mathbf{V}^{*}\|\right)\left\|\boldsymbol{\Xi}\right\|_{F}\\\
&\leq
2\left(\|\mathbf{V}_{\star}\widetilde{\mathbf{V}}_{\perp}^{*}\|+\|\mathbf{V}_{\star}\mathbf{V}_{\perp}^{*}\|\right)\left\|\boldsymbol{\Xi}\right\|_{F}\\\
&\leq
2\left(\frac{\|\mathsf{H}_{s}(\mathbf{X}^{(k)})-\mathbf{X}_{\star}\|}{\sigma_{r}(\mathbf{X}_{\star})-\sigma_{r+1}(\mathsf{H}_{s}(\mathbf{X}^{(k)}))}+\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{\sigma_{r}(\mathbf{X}_{\star})-\sigma_{r+1}(\mathbf{X}^{(k)})}\right)\left\|\boldsymbol{\Xi}\right\|_{F}\\\
&\leq
4\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{\sigma_{r}(\mathbf{X}_{\star})-\sigma_{r+1}(\mathbf{X}^{(k)})}\left\|\boldsymbol{\Xi}\right\|_{F}\leq
4\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|}{(1-1/48)\sigma_{r}(\mathbf{X}_{\star})}\left\|\boldsymbol{\Xi}\right\|_{F},\end{split}$

using again that
$\sigma_{r+1}(\mathbf{X}^{(k)})\leq\|\mathbf{X}_{\star}-\mathbf{X}^{(k)}\|\leq\frac{1}{48}\sigma_{r}(\mathbf{X}_{\star})$
due to Eq. 26.

Let now $\boldsymbol{\Xi}^{(k+1)}=\mathbf{X}^{(k+1)}-\mathbf{X}_{\star}$.
Combining Eq. 25 and Eq. 51-Eq. 52 we can proceed to estimate that

(55) $\displaystyle\begin{split}&\|\boldsymbol{\Xi}^{(k+1)}\|_{F}^{2}\\\ &\leq
c^{2}\left\|\mathbb{P}_{T_{k}^{\perp}}(\boldsymbol{\Xi}^{(k+1)})\right\|_{F}^{2}\\\
&\leq
2c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}\boldsymbol{\Xi}^{(k+1)}\big{\rangle}\\\
&+2\Big{[}\Big{(}\frac{8}{47/48}\Big{)}^{2}+\Big{(}\frac{4}{47/48}\Big{)}^{2}\Big{]}c^{2}\|\boldsymbol{\Xi}^{(k+1)}\|_{F}^{2}\frac{\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|^{2}}{\sigma_{r}^{2}(\mathbf{X}_{\star})}\\\
&\leq
2c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\boldsymbol{\Xi}^{(k+1)})\big{\rangle}+167c^{2}\|\boldsymbol{\Xi}^{(k+1)}\|_{F}^{2}\frac{1}{(19c)^{2}}\\\
&\leq
2c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\boldsymbol{\Xi}^{(k+1)})\big{\rangle}+\frac{1}{2}\|\boldsymbol{\Xi}^{(k+1)}\|_{F}^{2},\end{split}$

where the third inequality follows from Eq. 8 and Eq. 26. Hence, rearranging
LABEL:eq:B8intermediateStepII yields

$\displaystyle\|\boldsymbol{\Xi}^{(k+1)}\|_{F}^{2}\leq
4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\boldsymbol{\Xi}^{(k+1)})\big{\rangle}.$

By Lemma C.6, we know that $\mathbf{X}^{(k+1)}$ fulfills

$\displaystyle
0=\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k+1)})\big{\rangle}=\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\boldsymbol{\Xi}^{(k+1)})\big{\rangle}+\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{X}_{\star})\big{\rangle}$

such that we conclude that

$\displaystyle\left\|\boldsymbol{\Xi}^{(k+1)}\right\|^{2}\leq\left\|\boldsymbol{\Xi}^{(k+1)}\right\|_{F}^{2}$
$\displaystyle\leq
4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\boldsymbol{\Xi}^{(k+1)})\big{\rangle}$
$\displaystyle=-4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{X}_{\star})\big{\rangle}$
$\displaystyle=-4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\left(\big{\langle}\boldsymbol{\Xi}^{(k+1)},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\big{\rangle}+\big{\langle}\boldsymbol{\Xi}^{(k+1)},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\cdot\mathbf{X}_{\star}\big{\rangle}\right)$
$\displaystyle\leq
4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\\!\\!\left(\left\|W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\right\|_{*}\left\|\boldsymbol{\Xi}^{(k+1)}\right\|\\!\\!+\\!\left\|\boldsymbol{\Xi}^{(k+1)}\right\|\\!\left\|\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\cdot\mathbf{X}_{\star}\right\|_{1,2}\right)$
$\displaystyle=4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\left(\left\|W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\right\|_{*}\\!+\\!\left\|\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\cdot\mathbf{X}_{\star}\right\|_{1,2}\right)\left\|\boldsymbol{\Xi}^{(k+1)}\right\|,$

which completes the proof. We used in the penultimate line Hölder’s inequality
and that

$\displaystyle|\langle\mathbf{A},\mathbf{B}\rangle_{F}|=\left|\sum_{i,j}A_{i,j}B_{i,j}\right|\leq\sum_{i}\left\|\mathbf{A}_{i,:}\right\|_{2}\left\|\mathbf{B}_{i,:}\right\|_{2}\leq\Big{(}\max_{i}\left\|\mathbf{A}_{i,:}\right\|_{2}\Big{)}\cdot\sum_{i}\left\|\mathbf{B}_{i,:}\right\|_{2}\leq\left\|\mathbf{A}\right\|\left\|\mathbf{B}\right\|_{1,2},$

for all matrices $\mathbf{A},\mathbf{B}$.

###  D.5 Proof of Lemma C.10

Note that by Lemma C.2 and Eq. 28,
$S:=\mathrm{supp}(\mathsf{H}_{s}(\mathbf{X}^{(k)}))=\mathrm{supp}(\mathbf{X}_{\star})$.

Since by assumption $\delta_{k}\leq\rho_{s}(\mathbf{X}^{(k)})$ we have by
definition of $\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$ that
$\mathbf{Z}\leavevmode\nobreak\ :=\leavevmode\nobreak\
\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\leavevmode\nobreak\
\cdot\leavevmode\nobreak\ \mathbf{X}_{\star}$ is a matrix with row-support $S$
and rows

$\displaystyle\mathbf{Z}_{i,:}=\min\left\\{\frac{\delta_{k}^{2}}{\big{\|}(\mathbf{X}^{(k)})_{i,:}\big{\|}_{2}^{2}},1\right\\}(\mathbf{X}_{\star})_{i,:},=\frac{\delta_{k}^{2}}{\big{\|}(\mathbf{X}^{(k)})_{i,:}\big{\|}_{2}^{2}}(\mathbf{X}_{\star})_{i,:}$

for $i\in S$. Now note that if Eq. 28 holds, then

$\displaystyle\|\mathbf{Z}_{i,:}\|_{2}=\frac{\delta_{k}^{2}\|(\mathbf{X}_{\star})_{i,:}\|_{2}}{\|(\mathbf{X}^{(k)})_{i,:}\|_{2}^{2}}\leq\frac{\delta_{k}^{2}}{(1-\zeta)^{2}\rho_{s}(\mathbf{X}_{\star})},$

where we used in the last estimate that with Eq. 28 and
$\|(\mathbf{X}_{\star})_{i,:}\|_{2}\geq\rho_{s}(\mathbf{X}_{\star})$, for
$i\in S$, we have

$\displaystyle\|(\mathbf{X}^{(k)})_{i,:}\|_{2}^{2}$
$\displaystyle\geq(\|(\mathbf{X}_{\star})_{i,:}\|_{2}-\|(\mathbf{X}^{(k)})_{i,:}-(\mathbf{X}_{\star})_{i,:}\|_{2})^{2}\geq(\|(\mathbf{X}_{\star})_{i,:}\|_{2}-\zeta\rho_{s}(\mathbf{X}_{\star}))^{2}$
$\displaystyle\geq\Big{(}(1-\zeta)\rho_{s}(\mathbf{X}_{\star})\Big{)}\Big{(}(1-\zeta)\|(\mathbf{X}_{\star})_{i,:}\|_{2}\Big{)},$

for all $i\in S$. The claim easily follows since $\mathbf{Z}$ has only $s$
non-zero rows.

###  D.6 Proof of Lemma C.11

In the proof of Lemma C.11, we use a simple technical observation.

######  Lemma D.5.

Let $\mathbf{X}_{\star}\in\mathcal{M}_{r,s}$, let $\mathbf{X}^{(k)}$ be the
$k$-th iterate of Algorithm 1, and abbreviate
$\boldsymbol{\Xi}^{(k)}=\mathbf{X}^{(k)}-\mathbf{X}_{\star}$. Then the
following two statements hold true:

  1. 1.

If $\varepsilon_{k}\leq\sigma_{r+1}(\mathbf{X}^{(k)})$, then
$\varepsilon_{k}\leq\|\boldsymbol{\Xi}^{(k)}\|$.

  2. 2.

If $\delta_{k}\leq\rho_{s+1}(\mathbf{X}^{(k)})$, then
$\delta_{k}\leq\|\boldsymbol{\Xi}^{(k)}\|_{\infty,2}$.

######  Proof D.6.

By defining $[\mathbf{X}^{(k)}]_{r}$ to be the best rank-$r$ approximation of
$\mathbf{X}^{(k)}$ in any unitarily invariant norm, we bound

$\displaystyle\varepsilon_{k}\leq\sigma_{r+1}(\mathbf{X}^{(k)})=\|\mathbf{X}^{(k)}-[\mathbf{X}^{(k)}]_{r}\|\leq\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|=\|\boldsymbol{\Xi}^{(k)}\|,$

where the inequality follows the fact that $\mathbf{X}_{\star}$ is a rank-$r$
matrix.

Similarly, for the second statement, we have that

$\delta_{k}\leq\rho_{s+1}(\mathbf{X}^{(k)})=\|\mathbf{X}^{(k)}-\mathsf{H}_{s}(\mathbf{X}^{(k)})\|_{\infty,2}\leq\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|_{\infty,2}=\|\boldsymbol{\Xi}^{(k)}\|_{\infty,2},$

using that $\mathbf{X}_{\star}$ is $s$-row sparse.

######  Proof D.7 (of Lemma C.11).

First, we note that, using Lemma C.4, the observation in Remark C.8 yields

(56) $\begin{split}\|\boldsymbol{\Xi}\|_{F}^{2}&\leq
4c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\\\
&\qquad\cdot\Big{\langle}\boldsymbol{\Xi},\big{(}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}+\mathbb{P}_{S^{c}}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\big{)}\boldsymbol{\Xi}\Big{\rangle}\end{split}$

for all $\boldsymbol{\Xi}\in\mathbb{R}^{n_{1}\times n_{2}}$ as the assumption
Eq. 29 implies
$\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|\leq\min\Big{\\{}\frac{1}{48},\frac{1}{19c}\Big{\\}}\sigma_{r}(\mathbf{X}_{\star})$.
Thus, this holds in particular also for
$\boldsymbol{\Xi}^{(k)}=\mathbf{X}^{(k)}-\mathbf{X}_{\star}$. (Recall that
$\mathbf{U}$ and $\mathbf{V}$ contain the leading singular vectors of
$\mathbf{X}^{(k)}$, see Definition 2.1.) We estimate that

$\begin{split}&\sqrt{\big{\langle}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}^{(k)},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}\boldsymbol{\Xi}^{(k)}\big{\rangle}}=\left\|(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}})^{1/2}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\boldsymbol{\Xi}^{(k)}))\right\|_{F}\\\
&\leq\left\|(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}})^{1/2}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\mathbf{X}_{\star}))\right\|_{F}+\left\|(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}})^{1/2}(\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}(\mathbf{X}^{(k)}))\right\|_{F}\\\
&\leq\left\|(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}})^{1/2}(\mathbf{X}_{\star})\right\|_{F}+\sqrt{\sum_{i=r+1}^{\min(n_{1},n_{2})}\frac{\sigma_{i}^{2}}{\max\big{(}\frac{\sigma_{i}^{2}}{\varepsilon_{k}^{2}},1\big{)}}}\leq\left\|(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}})^{1/2}(\mathbf{X}_{\star})\right\|_{F}+\sqrt{n-r}\varepsilon_{k}\end{split}$

Furthermore, since
$\max(\varepsilon_{k},\|\boldsymbol{\Xi}^{(k)}\|)\leq\frac{1}{48}\sigma_{r}(\mathbf{X}_{\star})$
by assumption and $\varepsilon_{k}\leq\|\boldsymbol{\Xi}^{(k)}\|$ by Lemma
D.5, we can use a variant of Lemma C.9 to obtain

$\begin{split}\left\|(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}})^{1/2}(\mathbf{X}_{\star})\right\|_{F}&\leq\frac{48}{47}\left(\sqrt{r}\varepsilon_{k}+2\varepsilon_{k}\frac{\|\boldsymbol{\Xi}^{(k)}\|_{F}}{\sigma_{r}(\mathbf{X}_{\star})}+2\frac{\|\boldsymbol{\Xi}^{(k)}\|\|\boldsymbol{\Xi}^{(k)}\|_{F}}{\sigma_{r}(\mathbf{X}_{\star})}\right)\\\
&\leq
1.04\sqrt{r}\varepsilon_{k}+\frac{4.16\|\boldsymbol{\Xi}^{(k)}\|}{\sigma_{r}(\mathbf{X}_{\star})}\|\boldsymbol{\Xi}^{(k)}\|_{F}.\end{split}$

On the other hand, we note that $\boldsymbol{\Xi}^{(k)}$ restricted to $S^{c}$
coincides with the restriction of $\mathbf{X}^{(k)}$ to $S^{c}$ under
assumption Eq. 29, cf. Lemma C.2, and therefore

$\begin{split}\big{\langle}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}^{(k)},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\boldsymbol{\Xi}^{(k)}\big{\rangle}&=\big{\langle}\mathbb{P}_{S^{c}}\mathbf{X}^{(k)},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\mathbf{X}^{(k)}\big{\rangle}\\\
&=\sum_{i=s+1}^{n_{1}}\frac{\|(\mathbf{X}^{(k)})_{i,:}\|_{2}^{2}}{\max\\{\|(\mathbf{X}^{(k)})_{i,:}\|_{2}^{2}/\delta_{k}^{2},1\\}}\leq(n_{1}-s)\delta_{k}^{2}.\end{split}$

With the estimate of above, this implies that

$\displaystyle\begin{split}&\Big{\langle}\boldsymbol{\Xi},\big{(}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}+\mathbb{P}_{S^{c}}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\big{)}\boldsymbol{\Xi}\Big{\rangle}\\!\\\
&\leq\\!\Big{(}\\!\left\|(W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}})^{1/2}(\mathbf{X}_{\star})\right\|_{F}\\!\\!+\\!\sqrt{n-r}\varepsilon_{k}\Big{)}^{2}\\!\\!\\!+\\!(n_{1}-s)\delta_{k}^{2}\\\
&\leq\frac{13}{4}r\varepsilon_{k}^{2}+\frac{52\|\boldsymbol{\Xi}^{(k)}\|^{2}}{\sigma_{r}^{2}(\mathbf{X}_{\star})}\|\boldsymbol{\Xi}^{(k)}\|_{F}^{2}+3(n-r)\varepsilon_{k}^{2}+(n_{1}-s)\delta_{k}^{2}.\end{split}$

Inserting these estimates into Eq. 56, we obtain

$\|\boldsymbol{\Xi}^{(k)}\|_{F}^{2}\leq 4c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\left(\frac{13}{4}n\varepsilon_{k}^{2}+n_{1}\delta_{k}^{2}+\frac{52\|\boldsymbol{\Xi}^{(k)}\|^{2}}{\sigma_{r}^{2}(\mathbf{X}_{\star})}\|\boldsymbol{\Xi}^{(k)}\|_{F}^{2}\right).$

If now either one of the two equations
$\varepsilon_{k}=\sigma_{r+1}(\mathbf{X}^{(k)})$ or
$\delta_{k}=\rho_{s+1}(\mathbf{X}^{(k)})$ is true, it follows that

$\begin{split}\|\boldsymbol{\Xi}^{(k)}\|_{F}^{2}&\leq
c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\left(13n\varepsilon_{k}^{2}+4n_{1}\delta_{k}^{2}\right)+c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\frac{208\|\boldsymbol{\Xi}^{(k)}\|^{2}}{\sigma_{r}^{2}(\mathbf{X}_{\star})}\|\boldsymbol{\Xi}^{(k)}\|_{F}^{2}\\\
&\leq c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\left(13n\varepsilon_{k}^{2}+4n_{1}\delta_{k}^{2}\right)+\frac{1}{2}\|\boldsymbol{\Xi}^{(k)}\|_{F}^{2}\end{split}$

if the proximity condition
$\|\boldsymbol{\Xi}^{(k)}\|=\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|\leq\frac{1}{21c_{\left\|\mathcal{A}\right\|_{2\to
2}}}\sigma_{r}(\mathbf{X}_{\star})$ is satisfied. Rearranging the latter
inequality yields the conclusion of Lemma C.11.


# Continuous Real-Time Detection of Quasiparticle Trapping

J. T. Farmer  A. Zarassi  D. M. Hartsell  E. Vlachos  H. Zhang  E. M.
Levenson-Falk<EMAIL_ADDRESS>Department of Physics & Astronomy, Dornsife
College of Letters, Arts, & Sciences, University of Southern California

###### Abstract

Nonequilibrium quasiparticles are ubiquitous in superconducting electronics.
These quasiparticles can trap in the internal Andreev bound states of a phase-
biased Josephson junction, providing a mechanism for studying their presence
and behavior. We characterize a quasiparticle trapping detector device based
on a two-junction aluminum nanobridge superconducting quantum interference
device incorporated into a transmission-line resonator. When flux-biased,
distinct resonant frequencies develop depending on the trapped quasiparticle
number. We demonstrate continuous detection of up to 2 trapped quasiparticles
with signal to noise ratio of 1 in 60 ns. We discuss the possible optimization
and application of such detector devices.

††preprint: AIP/123-QED

Superconducting qubits and other low-temperature superconducting electronics
have ubiquitous populations of quasiparticles (QPs) far above their thermal
equilibrium prevalence Aumentado  _et al._ (2004); Martinis, Ansmann, and
Aumentado (2009); Catelani  _et al._ (2011); Rainis and Loss (2012); Wang  _et
al._ (2014); Vool  _et al._ (2014); Glazman and Catelani (2021). These QPs can
cause lossGustavsson  _et al._ (2016), spurious excitationSerniak  _et al._
(2018), and spectral noise when tunneling across qubit junctions, and can
induce correlated errors that are particularly difficult to address with error
correction algorithms Vepsäläinen  _et al._ (2020); McEwen  _et al._ . The QPs
may be generated by stray infrared photons Barends  _et al._ (2011), cosmic
rays or other high-energy radiation sources Vepsäläinen  _et al._ (2020);
Cardani  _et al._ (2021); McEwen  _et al._ ; Wilen  _et al._ (2021), or
materials defects Kurter  _et al._ . Many experiments have probed QP behavior
via their tunneling across Josephson junctions in charge-sensitive transmons
Sun  _et al._ (2012); Ristè  _et al._ (2013); Serniak  _et al._ (2018, 2019),
giving valuable insight into their effects on qubits. In these experiments
each tunneling measurement must be discrete, and they cannot distinguish
between 0 and 2 tunneling events. Trapping measurements Zgirski  _et al._
(2011); Levenson-Falk  _et al._ (2014); Hays  _et al._ (2020) provide a tool
for continuous, nonsaturating monitoring of QP behavior. Furthermore, QP traps
have been proposed as a way to mitigate QPs’ harmful effects on qubits Sun
_et al._ (2012); Wang  _et al._ (2014); Riwar  _et al._ (2016); Martinis
(2021), as QPs may diffuse great distances after being generated Patel  _et
al._ (2017); McEwen  _et al._ ; Wilen  _et al._ (2021). The trapping process
itself is thus worthy of study, in addition to providing insight into other QP
properties.

In this letter we characterize a device optimized for continuous, non-
saturating measurements of quasiparticle trapping in Andreev states. Using
standard microwave reflectometry we are able to distinguish 0, 1, and 2 or
more trapped quasiparticles in 60 ns with a signal-to-noise ratio (SNR) of 1.
We discuss straightforward improvements that can further improve SNR and allow
detection of many more trapped quasiparticles. Our device provides a prototype
for future detectors optimized for continuous measurements of quasiparticle
behavior and properties and for studies of the dynamics of Andreev states
coupled to resonant cavities.

In the semiconductor picture of the Josephson effect, supercurrent is carried
by quasiparticles/quasiholes traveling back and forth in 1 dimensional
conduction channels. At each boundary of the junction, the quasiparticle
(quasihole) reflects as a quasihole (quasiparticle) Andreev (1964). This
Andreev reflection process causes a $\pm 2e$ charge transfer, where $e$ is the
elementary charge, and thus transmits a Cooper pair across the junction. Each
conduction channel forms a pair of Andreev bound states with energies

$E_{A\pm}=\pm\Delta\sqrt{1-\tau\sin^{2}\frac{\delta}{2}}$ (1)

where $\Delta$ is the superconducting gap, $\tau$ is the tranmittivity of the
channel, and $\delta$ is the phase bias across the junction. The Fermi energy
is 0, so at temperatures $T\ll\Delta/k_{B}$, the positive-energy state is
normally unoccupied and the negative-energy state is normally occupied,
carrying the supercurrent across the junction. Quasiparticles in the bulk of
the superconductor can only exist at energies greater than $\Delta$, so such a
quasiparticle will find it energetically favorable to drop into the unoccupied
upper Andreev state, which carries the opposite supercurrent. This Andreev
trapping process “poisons” the channel, eliminating it from carrying
supercurrent, reducing the effective junction critical current, and increasing
the Josephson inductance.

Three dimensional aluminum nanobridge Josephson junctions achieve good phase
confinement and nonlinearity in an all-superconducting design Vijay  _et al._
(2009, 2010); Levenson-Falk, Vijay, and Siddiqi (2011). Importantly, these
junctions comprise a large number ($\sim 100-1000$) of conduction channels
with high transmittivities, approximately following the Dorokhov distribution
$\rho(\tau)=\frac{\pi\hbar G}{2e^{2}}\frac{1}{\tau\sqrt{1-\tau}}$, where $G$
is the junction’s normal-state conductance Dorokhov (1982). These junctions
thus form a large number of Andreev bound states that, when phase-biased,
function as quasiparticle traps. For an Al nanobridge with phase bias
$\delta\approx 0.47\pi$, a channel with $\tau\approx 1$ has $E_{A+}/h\approx
30.4$ GHz, or a trap depth of $(\Delta-E_{A})/h\approx 10.7$ GHz, far greater
than the thermal energy at 15 mK. These junctions thus function as
quasiparticle traps with a large number of deep trap states when they are
phase-biased.

Our device consists of a co-planar waveguide (CPW) resonator in which the
center trace is terminated by a pair of Aluminum nanobridge junctions forming
a superconducting quantum interference device (SQUID), as shown in Fig. 1
panels (a) and (b). To fabricate the device, we use photolithography and
electron-beam evaporation to deposit the open ended CPW trace on a high
resistivity Si substrate. We then use electron-beam lithography and another
round of evaporation to deposit the junctions with connection pads overlaying
the CPW trace and ground plane. Finally, photolithography is used to open
windows in the resist for ion milling and deposition of superconducting
contacts between the SQUID and CPW. Fabrication details are further explained
in the supporting informationsup . The fundamental (quarter-wavelength) mode
of our resonator with zero flux through the SQUID is at
$\omega_{0}(\phi=0)=2\pi\times 4.302$ GHz and has linewidth $\kappa=2\pi\times
250$ kHz, largely set by the coupling to the microwave feedline.

Figure 1: [color online] (a) Device schematic. A CPW resonator (green) is
grounded via a SQUID (magenta) formed by two superconducting Al nanobridge
junctions. Flux bias through the SQUID phase-biases the junctions. (b) Optical
image of the CPW resonator with inset SEM images of the SQUID (magenta) and a
nanobridge junction (orange). (c) Simplified measurement schematic. A tone at
$\omega_{d}$ continuously drives the flux tunable resonator. All input and
output signals are filtered at base stage by K&L 12 GHz and custom Eccosorb
110 low-pass filters. The reflected signal is amplified by a TWPA at base
stage followed by a HEMT at 4K and room temperature amplifiers. Note that
isolators between HEMT and room temperature amplifiers are not shown to
conserve space. The amplified signal is homodyne demodulated and the in phase
(I) and quadrature (Q) components are digitized by an Alazar ATS9371 at 300
MSa/s after 15 MHz low-pass filtering. The recorded signal is then down-
sampled to 10 MSa/s before saving.

Our measurement setup is shown schematically in Fig. 1 (c). A signal generator
provides a drive tone at $\omega_{d}(\phi)\approx\omega_{0}(\phi)-\kappa/2$. A
power splitter sends half of this power into the dilution refrigerator where
it is attenuated then filtered by K&L 12 GHz and Eccosorb low pass filters.
The drive tone is circulated to reflect off our device, which is flux tunable
via a DC coil in the packaging, and amplified by a travelling wave parametric
amplifierMacklin  _et al._ (2015) (TWPA) which is pumped at 8.078 GHz. The
reflected signal is further amplified before IQ demodulation with reference to
the original signal. The in-phase (I) and quadrature (Q) components are 15 MHz
low pass filtered before digitization at 300 MSa/s. Raw data is down-sampled
to 10 MSa/s before saving.

We first characterize our device with ensemble-averaged VNA measurements of
the resonance. Fig. 2 (a) shows resonance measurements at flux biases of
$\phi=0$ and $\phi=0.49$, taken on a previous cooldown in which the K&L and
Eccosorb filters were removed from the setup. In this configuration QP
trapping is more frequent, which we attribute to more infrared radiation
reaching the device via the unfiltered lines. The $\phi=0.49$ trace in orange
shows two shallow peaks at $\sim$ 0.5 and 1 MHz below the main resonance.
These are the resonance peaks with 1 and 2 trapped QPs, respectively, showing
the resonance shifting due to the change in nanobridge inductance. These
ensemble measurements average over all possible QP trapping configurations,
and so a resonance peak amplitude corresponds to the probability of that
configuration. We then move on to time-domain IQ measurements as described
above. Panels (b) and (c) are log-scale histograms showing 30 seconds of
continuous IQ data for $\phi=0$ (in blue) and $\phi=0.47$ (in orange),
respectively. This data has been convolved with a Gaussian window of effective
integration time $\int_{-\infty}^{\infty}f(t)dt=3~{}\mu$s. In the finite-flux
data of panel (c) we can immediately see 3 distinct modes with excellent
separation in the log scale plots. The darkest peak (lower left), with the
most counts by far, is due to the resonance when we have 0 QPs in Andreev
traps. The next darkest (upper center) is from having a single QP trapped
while the lightest mode (middle right) is from 2 or more QPs in Andreev traps.
Trapping of more than 2 quasiparticles moves the resonance multiple linewidths
and thus saturates the change in response to additional trapping. We later
discuss methods for avoiding this saturation. We have verified that these
modes are indeed due to Andreev trapping of quasiparticles by measuring the
weights of each mode in the presence of a “clearing tone” at 17 GHz sup . We
find that the 1- and 2-quasiparticle mode counts decrease while the
0-quasiparticle mode counts increase. We have also observed that the
separations and weights of modes both increase as a function of flux, which
agrees qualitatively with a quasiparticle trapping picture sup .

Figure 2: [color online] (a) Ensemble measurement of resonator response at
nearly half a flux quantum (orange) and zero flux (blue). Note that for
$\phi=0.49$, distinct peaks are visible near -0.5 and -1 MHz below resonance,
corresponding to 1 and 2 trapped QPs, respectively. (b) Histogram of
continuous IQ data taken at 0 flux for 10 s. Data has 10 MHz sample rate and
has been convolved with a gaussian window with effective integration time of 3
$\mu$s. (c) Data taken at $\phi=0.47$ with the same procedure as panel b. The
darkest mode is due to the resonance with 0 trapped QPs. The second darkest,
located near (I,Q) = (12 mV, 15 mV), corresponds to 1 trapped QP and the last
mode corresponds to 2 or more trapped QPs as this mode corresponds to the
resonance moving far from the drive frequency.

We now turn to extracting the QP trap occupation as a function of time from
the continuous IQ data. We first use a Gaussian mixture expectation-
maximization algorithm implemented by the python module available from scikit-
learnPedregosa  _et al._ (2011). This module takes in a subset of data,
assigns each point to one of the specified number of modes, then tweaks
assignments and mode parameters to maximize the total likelihood for all data
and all modes. The result is a set of 3 Gaussian modes describing the data.
These modes are shown by their 1-$\sigma$ (solid) and 2-$\sigma$ (dashed)
contours for each mode overlying the histogram in Fig. 3 (a). Each time-series
point is then assigned to the mode with the highest posterior probability
after Bayesian updating with a 3 sample rolling memory windowsup . However, we
find the Gaussian mixture module to be unreliable in terms of quality of fit
and reproducibility of Gaussian mode parameters, varying significantly even
when refitting the same data with the same initial guess. This is immediately
apparent from the assigned distributions shown, which do not faithfully
represent the means of the 1- and 2-QP modes.

Figure 3: [color online] (a) Initial clustering of data using the scikit-learn
Gaussian mixture module produces modes with 1$\sigma$ (solid) and 2$\sigma$
(dashed) contours for 0, 1, and 2 or more trapped QPs in light blue, dark
blue, and orange, respectively. (b) To improve extraction fidelity, subsets of
the data in which the occupation stays in a given mode for a substantial
amount of time ($4/\Gamma_{i}$) are individually fit to a Gaussian
distribution. The means and covariances of each mode are then fixed and a
final fit to the full dataset is performed with the only free parameters being
the weight of each mode. (c) A segment of data plotted as a time series with I
in brown and Q in dashed magenta. The trap occupation, extracted as the
maximum posterior probability associated with the modes of panel b for each
point in data (after a 3 sample window of Bayesian updating), is encoded in
the background color as light blue, dark blue, and orange for 0, 1, and 2+
QPs, respectively.

To improve the reliability of the fit, we need to “initialize” the trapped
quasiparticle configuration, thereby isolating a single Gaussian mode for
independent fitting. This is a challenge, as we have no direct control of the
actual trapped quasiparticle number. Fortunately, the Gaussian mixture module
does a good job of assigning the majority of each mode to the correct
occupation. We use this initial assignment to fit the mean decay rate of each
mode $\Gamma_{i}$ and then identify periods in the time series when the
extracted QP occupation is stationary for a substantial period of time, at
least $4/\Gamma_{i}$. By stitching many of these “quiet periods” together, we
build up large distributions of data points that are all pre-assigned to a
mode. These “quiet period” distributions are then independently fit to
Gaussians to fix their means and covariances. Finally, the full dataset is fit
to a mixture of Gaussians with these same means and covariances, with the only
free parameter being the weight of each mode. Fig. 3 (b) shows the result of
this procedure, with 1-$\sigma$ (solid) and 2-$\sigma$ (dashed) contours
overlaying the data histogram after this procedure. It is immediately evident
that our “quiet periods” method produces a better quality of fit. Extraction
of the QP occupation is performed identically to previous method with updated
Gaussian parameters. Using this new assignment algorithm we proceed to assign
each time point to a QP occupation. Fig. 3 (c) shows 1 ms of data with the
same 3 $\mu$s Gaussian convolution. The background shading represents the
extracted QP trap occupation (light blue for 0 QPs, dark blue for 1, and
orange for 2 or more). Clear transitions between the 3 configurations are
immediately visible and appear to be faithfully captured by the assignment
algorithm. We find our signal-to-noise ratio at our operating parameters
($\sim$ 16 photon drive power, 3 $\mu$s integration), defined here as the
separation of mode centers in the IQ plane squared divided by product of
covariances of each mode projected onto the line connecting them, to be 49 for
$0-1$ distinguishability, 49 for $1-2$, and 44 for $0-2$. Extrapolating from
the $0-1$ SNR, we find that our setup has a noise floor of $3.5\times 10^{-4}$
QPs /$\sqrt{\text{Hz}}$, i.e., we can detect 0.00035 of the resonance shift
from a QP trapping with SNR = 1 after a half second of integration, assuming a
stationary occupation.

We note that our device is not fully optimized for high sensitivity. While the
resonant frequency was kept low in order to avoid interactions with the
Andreev gap frequencies, raising the resonance slightly to $\sim$ 6 GHz by
shortening the waveguide would increase the participation ratio of the
nanobridge inductance to total inductance while still remaining far below the
gap energy. Similarly, reducing the resonator’s characteristic impedance from
50 $\Omega$ to an easily-achievable $\sim 25\Omega$ would further reduce the
linear inductance, increasing sensitivity. We also note that the TWPA used as
a first-stage amplifier had a moderate $\approx$ 15 dB gain (due to being
operated near the edge of its bandwidth) and does not have a quantum-limited
noise temperature. Simply adding a standard parametric amplifier with a near-
quantum-limited noise temperature and 20 dB of gain as a preamplifier should
further improve our SNR. Future work will take advantage of these SNR
improvements. We also note the possibility of detecting higher numbers of
trapped QPs by probing the device with multiple probe tones simultaneously. By
performing heterodyne measurement of probe tones centered on, for example, the
1-, 3-, and 5-QP resonant frequencies and correlating the measured outcomes,
we should be able to detect up to 6 trapped QPs with a similar device.

In conclusion, we have developed a device for the ultra-low-noise continuous
detection of up to 2 quasiparticles trapping in Andreev bound states. Our
device is capable of detecting a trapped QP with SNR of 1 in 60 ns, giving it
a noise floor of $3.5\times 10^{-4}$ QPs /$\sqrt{\text{Hz}}$. Straightforward
extensions are possible to lower noise floors and detection of higher numbers
of trapped QPs. Our device can be used for QP studies including statistical
analysis of trapping and untrapping rates as well as trap occupation,
spectroscopic measurements of the trapped QP energy distributions, and testing
of QP mitigation techniques.

The authors would like to acknowledge J. Aumentado and L. Glazman for useful
discussions. We acknowledge MIT Lincoln Laboratory and IARPA for providing the
TWPA used in this work. This work was funded by the AFOSR YIP under grant
FA9550-19-1-0060 and by the NSF DMR under grant DMR-1900135. The data that
support the findings of this study are available upon request.

## References

  * Aumentado  _et al._ (2004) J. Aumentado, M. W. Keller, J. M. Martinis, and M. H. Devoret, “Nonequilibrium Quasiparticles and 2e Periodicity in Single-Cooper-Pair Transistors,” Phys. Rev. Lett. 92, 066802 (2004).
  * Martinis, Ansmann, and Aumentado (2009) J. M. Martinis, M. Ansmann, and J. Aumentado, “Energy decay in superconducting josephson-junction qubits from nonequilibrium quasiparticle excitations,” Phys. Rev. Lett. 103, 097002 (2009).
  * Catelani  _et al._ (2011) G. Catelani, J. Koch, L. Frunzio, R. J. Schoelkopf, M. H. Devoret, and L. I. Glazman, “Quasiparticle relaxation of superconducting qubits in the presence of flux,” Phys. Rev. Lett. 106, 077002 (2011).
  * Rainis and Loss (2012) D. Rainis and D. Loss, “Majorana qubit decoherence by quasiparticle poisoning,” Phys. Rev. B. 85, 174533 (2012).
  * Wang  _et al._ (2014) C. Wang, Y. Y. Gao, I. M. Pop, U. Vool, C. Axline, T. Brecht, R. W. Heeres, L. Frunzio, M. H. Devoret, G. Catelani, L. I. Glazman, and R. J. Schoelkopf, “Measurement and control of quasiparticle dynamics in a superconducting qubit,” Nat. Commun. 5 (2014), 10.1038/ncomms6836.
  * Vool  _et al._ (2014) U. Vool, I. M. Pop, K. Sliwa, B. Abdo, C. Wang, T. Brecht, Y. Y. Gao, S. Shankar, M. Hatridge, G. Catelani, M. Mirrahimi, L. Frunzio, R. J. Schoelkopf, L. I. Glazman, and M. H. Devoret, “Non-poissonian quantum jumps of a fluxonium qubit due to quasiparticle excitations,” Phys. Rev. Lett. 113 (2014), 10.1103/PhysRevLett.113.247001.
  * Glazman and Catelani (2021) L. Glazman and G. Catelani, “Bogoliubov quasiparticles in superconducting qubits,” SciPost Phys. Lect. Notes 31, 031 (2021).
  * Gustavsson  _et al._ (2016) S. Gustavsson, F. Yan, G. Catelani, J. Bylander, A. Kamal, J. Birenbaum, D. Hover, D. Rosenberg, G. Samach, A. P. Sears, S. J. Weber, J. L. Yoder, J. Clarke, A. J. Kerman, F. Yoshihara, Y. Nakamura, T. P. Orlando, and W. D. Oliver, “Suppressing relaxation in superconducting qubits by quasiparticle pumping,” Science 354, 1573–1577 (2016).
  * Serniak  _et al._ (2018) K. Serniak, M. Hays, G. De Lange, S. Diamond, S. Shankar, L. D. Burkhart, L. Frunzio, M. Houzet, and M. H. Devoret, “Hot nonequilibrium quasiparticles in transmon qubits,” Phys. Rev. Lett. 121, 157701 (2018).
  * Vepsäläinen  _et al._ (2020) A. P. Vepsäläinen, A. H. Karamlou, J. L. Orrell, A. S. Dogra, B. Loer, F. Vasconcelos, D. K. Kim, A. J. Melville, B. M. Niedzielski, J. L. Yoder, S. Gustavsson, J. A. Formaggio, B. A. VanDevender, and W. D. Oliver, “Impact of ionizing radiation on superconducting qubit coherence,” Nature 584, 551–556 (2020).
  * (11) M. McEwen, L. Faoro, K. Arya, A. Dunsworth, T. Huang, S. Kim, B. Burkett, A. Fowler, F. Arute, J. C. Bardin, A. Bengtsson, A. Bilmes, B. B. Buckley, N. Bushnell, Z. Chen, R. Collins, S. Demura, A. R. Derk, C. Erickson, M. Giustina, S. D. Harrington, S. Hong, E. Jeffrey, J. Kelly, P. V. Klimov, F. Kostritsa, P. Laptev, A. Locharla, X. Mi, K. C. Miao, S. Montazeri, J. Mutus, O. Naaman, M. Neeley, C. Neill, A. Opremcak, C. Quintana, N. Redd, P. Roushan, D. Sank, K. J. Satzinger, V. Shvarts, T. White, Z. J. Yao, P. Yeh, J. Yoo, Y. Chen, V. Smelyanskiy, J. M. Martinis, H. Neven, A. Megrant, L. Ioffe, and R. Barends, “Resolving catastrophic error bursts from cosmic rays in large arrays of superconducting qubits,”  arXiv:2104.05219 (2021) .
  * Barends  _et al._ (2011) R. Barends, J. Wenner, M. Lenander, Y. Chen, R. C. Bialczak, J. Kelly, E. Lucero, P. O’Malley, M. Mariantoni, D. Sank, H. Wang, T. C. White, Y. Yin, J. Zhao, A. N. Cleland, J. M. Martinis, and J. J. Baselmans, “Minimizing quasiparticle generation from stray infrared light in superconducting quantum circuits,” Applied Physics Letters 99, 113507 (2011).
  * Cardani  _et al._ (2021) L. Cardani, F. Valenti, N. Casali, G. Catelani, T. Charpentier, M. Clemenza, I. Colantoni, A. Cruciani, G. D’Imperio, L. Gironi, L. Grünhaupt, D. Gusenkova, F. Henriques, M. Lagoin, M. Martinez, G. Pettinari, C. Rusconi, O. Sander, C. Tomei, A. V. Ustinov, M. Weber, W. Wernsdorfer, M. Vignati, S. Pirro, and I. M. Pop, “Reducing the impact of radioactivity on quantum circuits in a deep-underground facility,” Nat. Commun. 12 (2021).
  * Wilen  _et al._ (2021) C. D. Wilen, S. Abdullah, N. A. Kurinsky, C. Stanford, L. Cardani, G. D’Imperio, C. Tomei, L. Faoro, L. B. Ioffe, C. H. Liu, A. Opremcak, B. G. Christensen, J. L. DuBois, and R. McDermott, “Correlated charge noise and relaxation errors in superconducting qubits,” Nature 594, 369–373 (2021).
  * (15) C. Kurter, C. E. Murray, R. T. Gordon, B. B. Wymore, M. Sandberg, R. M. Shelby, A. Eddins, V. P. Adiga, A. D. K. Finck, E. Rivera, A. A. Stabile, B. Trimm, B. Wacaser, K. Balakrishnan, A. Pyzyna, J. Sleight, M. Steffen, and K. Rodbell, “Quasiparticle tunneling as a probe of Josephson junction quality and capacitor material in superconducting qubits,”  arXiv:2106.11488 (2021) .
  * Sun  _et al._ (2012) L. Sun, L. Dicarlo, M. D. Reed, G. Catelani, L. S. Bishop, D. I. Schuster, B. R. Johnson, G. A. Yang, L. Frunzio, L. Glazman, M. H. Devoret, and R. J. Schoelkopf, “Measurements of quasiparticle tunneling dynamics in a band-gap-engineered transmon qubit,” Phys. Rev. Lett. 108, 230509 (2012).
  * Ristè  _et al._ (2013) D. Ristè, C. C. Bultink, M. J. Tiggelman, R. N. Schouten, K. W. Lehnert, and L. Dicarlo, “Millisecond charge-parity fluctuations and induced decoherence in a superconducting transmon qubit,” Nat. Commun. 4, 1–6 (2013).
  * Serniak  _et al._ (2019) K. Serniak, S. Diamond, M. Hays, V. Fatemi, S. Shankar, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, “Direct Dispersive Monitoring of Charge Parity in Offset-Charge-Sensitive Transmons,” Phys. Rev. Applied 12, 014052 (2019).
  * Zgirski  _et al._ (2011) M. Zgirski, L. Bretheau, Q. Le Masne, H. Pothier, D. Esteve, and C. Urbina, “Evidence for long-lived quasiparticles trapped in superconducting point contacts,” Phys. Rev. Lett. 106, 257003 (2011).
  * Levenson-Falk  _et al._ (2014) E. M. Levenson-Falk, F. Kos, R. Vijay, L. Glazman, and I. Siddiqi, “Single-quasiparticle trapping in aluminum nanobridge josephson junctions,” Phys. Rev. Lett. 112, 047002 (2014).
  * Hays  _et al._ (2020) M. Hays, V. Fatemi, K. Serniak, D. Bouman, S. Diamond, G. de Lange, P. Krogstrup, J. Nygård, A. Geresdi, and M. H. Devoret, “Continuous monitoring of a trapped superconducting spin,” Nat. Phys. 16, 1103–1107 (2020).
  * Riwar  _et al._ (2016) R.-P. Riwar, A. Hosseinkhani, L. D. Burkhart, Y. Y. Gao, R. J. Schoelkopf, L. I. Glazman, and G. Catelani, “Normal-metal quasiparticle traps for superconducting qubits,” Phys. Rev. B 94, 104516 (2016).
  * Martinis (2021) J. M. Martinis, “Saving superconducting quantum processors from decay and correlated errors generated by gamma and cosmic rays,” npj Quantum Information 7, 1–9 (2021).
  * Patel  _et al._ (2017) U. Patel, I. V. Pechenezhskiy, B. L. T. Plourde, M. G. Vavilov, and R. McDermott, “Phonon-mediated quasiparticle poisoning of superconducting microwave resonators,” Phys. Rev. B 96, 220501 (2017).
  * Andreev (1964) A. F. Andreev, “Thermal conductivity of the intermediate state in superconductors,” JETP 19, 1228 (1964).
  * Vijay  _et al._ (2009) R. Vijay, J. D. Sau, M. L. Cohen, and I. Siddiqi, “Optimizing anharmonicity in nanoscale weak link Josephson junction oscillators,” Phys. Rev. Lett. 103, 087003 (2009).
  * Vijay  _et al._ (2010) R. Vijay, E. M. Levenson-Falk, D. H. Slichter, and I. Siddiqi, “Approaching ideal weak link behavior with three dimensional aluminum nanobridges,” Applied Physics Letters 96, 3 (2010).
  * Levenson-Falk, Vijay, and Siddiqi (2011) E. M. Levenson-Falk, R. Vijay, and I. Siddiqi, “Nonlinear microwave response of aluminum weak-link Josephson oscillators,” Applied Physics Letters 98, 3 (2011).
  * Dorokhov (1982) O. Dorokhov, “Transmission coefficient and the localization length of an electron in N bound disordered chains,” Jetp Lett. 36, 259 (1982).
  * (30) See supplementary material.
  * Macklin  _et al._ (2015) C. Macklin, K. O’brien, D. Hover, M. Schwartz, V. Bolkhovsky, X. Zhang, W. Oliver, and I. Siddiqi, “A near–quantum-limited josephson traveling-wave parametric amplifier,” Science 350, 307–310 (2015).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” Journal of Machine Learning Research 12, 2825–2830 (2011).
  * Kulik and Omel’Yanchuk (1975) I. Kulik and A. Omel’Yanchuk, “Contribution to the microscopic theory of the josephson effect in superconducting bridges,” Soviet Journal of Experimental and Theoretical Physics Letters 21, 96 (1975).

##  I Supporting Information to “Continuous Real-Time Detection of
Quasiparticle Trapping”

###  I.1 Fabrication details

We start with a high-resistivity intrinsic silicon wafer with 5/50 nm Ti/Au
markers bordering a 1 cm square chip. AZ1512 photoresist is spun on at 4000
RPM for 60 seconds, followed by a pre-exposure bake at 100∘C for 60 seconds.
We use shadow-mask photolithography to pattern the coplanar waveguide
resonator and the chip is developed in AZ400K at room temperature for 45
seconds. Aluminum is deposited via e-beam evaporation to a thickness of 100
nm. Liftoff is performed by soaking in acetone for 30 minutes followed by a 30
second sonication.

To fabricate the junctions, we next spin the chip with MMA EL6 at 3200 RPM for
60 seconds followed by a 5 minute bake at 170∘C. This results in a layer of
MMA around 100 nm thick. We then spin a layer of PMMA A2 at 3600 RPM followed
by a 5 minute bake at 180∘C to achieve a 60 nm layer of PMMA on top. The SQUID
is patterned with e-beam lithography using a beam current of $\sim$18 pA and a
2 nm step size, with some intentionally low dose elements on one side of the
pattern that define undercuts – regions where the MMA develops away but a lip
of PMMA remains. The junctions are defined as 100 nm long lines bridging the
two sides of the SQUID. The chip is developed in a -15∘C chilled beaker of
MIBK:IPA 1:3 for 60 seconds, followed immediately by drying under a strong jet
of N2. The timing is crucial as any delay between development and drying will
warm the developer and result in overdeveloped features. Aluminum is deposited
via e-beam evaporation to 8 nm at direct incidence, then the substrate stage
is tilted to 35 degrees and another 62 nm thickness is accumulated on chip.
This angled evaporation is necessary to define the 3 dimensional nanobridge as
the 35 degree incidence blocks the line of sight to the substrate in the
junction line while allowing it in the SQUID leads. Liftoff is performed by a
2 hour soak in acetone warmed to 45∘C, followed by a 10 second sonication.

The last step is to make superconducting contact between the resonator and the
SQUID. AZ1512 is spun on as before and photolithography defines a pair of
rectangles for each device. This is developed in AZ400K as before. Each pair
of rectangles opens a window that covers a portion of the SQUID contact pad
and the resonator. The chip is then loaded into the e-beam evaporator and an
end-Hall ion source running Argon is used to remove surface oxides from the
SQUID and resonator portions exposed in these windows. Our recipe is 80 V at
3.5 A discharge (resulting in $\sim$50 eV neutral Argon bombardment) at 35
degree incidence for 270 seconds. We find that this is enough to remove the
native oxide layer while not etching through the thinnest (70 nm) layer of
Aluminum. Without breaking vacuum, we deposit 250 nm of Al to ensure good
contact between resonator and SQUID. Liftoff is performed by a 30 minute soak
in acetone followed by 10 second sonication.

###  I.2 Fitting flux tuning and participation ratio

To characterize the junction properties we fit the device resonant frequency,
$\omega_{0}$, as a function of applied flux, $\phi$, and extract the Josephson
inductance participation ratio at zero flux, $q_{0}$. To fit this, we need to
derive the flux tuning equation. We assume we have two symmetric junctions
which each approximately obey the KO-1 current-phase relation (CPR) of an
ideal short diffusive metallic weak link Kulik and Omel’Yanchuk (1975),

$I(\delta)=\frac{\Delta
N_{e}}{4\varphi_{0}}\cos\frac{\delta}{2}\tanh^{-1}(\sin\frac{\delta}{2}),$
(S1)

where $\Delta$ is the superconducting gap, $N_{e}$ is the effective number of
channels, $\varphi_{0}$ is the reduced flux quantum, and $\delta=\pi\phi$ is
the phase bias across the junction. This gives a SQUID inductance,

$L_{S}(\phi)=\frac{L_{j,0}}{2}[1-\sin\frac{\pi\phi}{2}\tanh^{-1}(\sin\frac{\pi\phi}{2})],$
(S2)

where the Josephson inductance per bridge at zero flux is

$L_{j,0}\equiv\frac{8\varphi_{0}^{2}}{\Delta N_{e}}.$ (S3)

With the definition of the Josephson inductance participation ratio as

$q(\phi)=\frac{L_{S}(\phi)}{L+L_{S}(\phi)}$ (S4)

and $q_{0}\equiv q(0)$, we can rewrite the resonant frequency function as

$\omega(\phi)=\omega_{0}\left[1+q_{0}\frac{\sin\frac{\pi\phi}{2}\tanh^{-1}(\sin\frac{\pi\phi}{2})}{1-\sin\frac{\pi\phi}{2}\tanh^{-1}(\sin\frac{\pi\phi}{2})}\right]^{-1/2}.$
(S5)

Figure S1: The extracted resonant frequency as a function of flux. The solid
curve is the fit to the equation in the supplemental text. We extract a
Josephson inductance participation ratio at zero flux of 0.01341.

We measure $S_{11}$ as a function of applied flux to extract the resonant
frequency, then fit it to the above equation. The results for the device in
the main text are presented in Fig. S1. We find our Josephson participation
ratio at zero flux is $q_{0}\approx 0.01341$. We obtain the resonator’s linear
inductance of $L=1.83$ nH via finite element EM simulations. From this and the
participation ratio we extract the SQUID inductance at zero flux as
$L_{S}(0)\approx 24$ pH, corresponding to an effective channel count
$N_{e}\approx 660$ per junction.

We note that our 100 nm long nanobridges have CPRs which are not quite equal
to the KO-1 formula Vijay  _et al._ (2009, 2010); Levenson-Falk, Vijay, and
Siddiqi (2011). However, we have found these junctions’ CPRs are very well
approximated by an ideal KO-1 junction in series with a small linear
inductance. This approximation leads to an insignificant modification of the
participation ratio, as this small extra linear inductance simply adds to the
large resonator linear inductance, and so we neglect it.

###  I.3 Quiet periods Gaussian fit and occupation extraction

We first fit the IQ voltage-vs-time data using the scikit-learn Gaussian
Mixture module described in the main text to extract a rough occupation. The
next task is to separate out subsets of data in which the mode (i. e. the
trapped QP number) is stationary for a significant portion of time. To do
this, we first fit the lifetimes of each occupation by building a histogram of
the times spent in each mode. The portion of this histogram above the mean is
fit to an exponential to extract the lifetime of the mode. For each mode, we
separate out only those subsets of data in which the mode is constant for
longer than 4 times that mode’s lifetime. We trim the first and last 5% of the
lifetime from each subset before stitching them together to exclude any
ringdown effects. This provides 3 datasets, one for each mode, where we can
confidently assign the data to a certain QP occupation.

Each of these datasets is histogrammed and independently fit to a 2D Gaussian
distribution. Using the means and covariances of each fit as fixed values, we
perform another fit of a histogram of the original full dataset to a sum of
all 3 Gaussians with their relative weights as the only free parameters. This
final step is of critical importance as the weights of each mode are a
significant contributor ($P(M_{i})$) to the calculation of posterior
probabilities used in extracting the occupation. Since we expect the resonator
to remain in the same mode for many consecutive points in the time series, we
apply a Bayesian inference scheme with a memory window of 3 points. i.e., the
probability that the $k_{th}$ point belongs to the $i_{th}$ mode is

$\displaystyle
P(M_{i}|\vec{x}_{k},\vec{x}_{k-1},\vec{x}_{k-2})=\frac{P(\vec{x}_{k}|M_{i})P(M_{i}|\vec{x}_{k-1},\vec{x}_{k-2})}{P(\vec{x}_{k}|M_{i})P(M_{i}|\vec{x}_{k-1},\vec{x}_{k-2})+\neg
P(\vec{x}_{k}|M_{i})\neg P(M_{i}|\vec{x}_{k-1},\vec{x}_{k-2})}$

where

$\displaystyle
P(M_{i}|\vec{x}_{k-1},\vec{x}_{k-2})=\frac{P(\vec{x}_{k-1}|M_{i})P(M_{i}|\vec{x}_{k-2})}{P(\vec{x}_{k}|M_{i})P(M_{i}|\vec{x}_{k-2})+\neg
P(\vec{x}_{k}|M_{i})\neg P(M_{i}|\vec{x}_{k-2})},$

and $P(M_{i}|\vec{x}_{k-2})$ is the posterior probability of the Gaussian mode
$M_{i}$ conditioned on the data, defined as

$P(M_{i}|\vec{x}_{k})=\frac{P(\vec{x}_{k}|M_{i})P(M_{i})}{\sum_{j=0}^{2}P(\vec{x}_{k}|M_{j})P(M_{j})},$
(S6)

Finally, the quasiparticle occupation at the $k_{th}$ point in the time series
is taken as the mode with maximum
$P(M_{i}|\vec{x}_{k},\vec{x}_{k-1},\vec{x}_{k-2})\;\forall\;i\in\\{0,1,2\\}.$

###  I.4 Clearing tone

As described in the main text, we expect the trap depth to be $\Delta-
E_{A}(\phi=0.47)\approx 10.7$ GHz. To test that our observed behavior is due
to quasiparticles trapping in Andreev states, we attempt to clear these
trapped quasiparticles by injecting a tone with sufficient energy to excite
the trapped quasiparticles into the continuum of available states above the
gap. Fig. S2 shows the microwave response of the resonator as a function of a
17 GHz tone drive power. Ensemble measurements of the reflection coefficient
show that, when driven strongly, the resonance curve recovers the single
Lorentzian shape that it has at 0 flux. We also see that the second and third
modes apparent in the continuous IQ measurements have greatly reduced
amplitudes, indicating that QPs spend far less time in trap states as the 17
GHz tone clears them out. Combined with the qualitative flux dependence of
trapping —trapping modes move farther apart and become stronger as the flux
bias becomes deeper— we take this as definitive evidence of QP trapping,
similar to that shown in ensemble measurements of similar devices Levenson-
Falk  _et al._ (2014).

Figure S2: Microwave response of the resonance at flux bias 0.45 as a function
of clearing tone power. The legend in (a) shows the approximate power of the
clearing tone at the plane of the device. The red curve is at low power and
shows a significant bump on the low frequency side due to averaging over many
configurations of quasiparticle occupation. At high clearing tone power (blue
curve) we see the resonator response narrows and appears to be more over-
coupled. We stress that the ratio of internal loss to external coupling is not
actually changing, rather, the response is taller because we are no longer
averaging over other configurations of quasiparticle occupation. (b) Histogram
of time series data with the clearing tone at the same power as the red trace
in (a). (c) Histogram of time series data with the clearing tone at the same
power as the blue trace in (a). Note that any trapped quasiparticles are
quickly excited back above the gap, greatly reducing the weights of the 1- and
2-QP modes and increasing the weight of the 0-QP mode.